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ABSTRACT 

This  work  is  aimed  at  developing  a  better  theoretical  under¬ 
standing  of  the  principal  transport  processes,  namely  diffusion  and 
fluid  flow,  when  occurring  within  unconsolidated  porous  media,  in 
particular  when  occurring  within  an  homogeneous  swarm  of  spheres. 

The  principal  objective  has  been  to  develop  an  efficient  geometric 
model  to  represent  such  a  system,  a  model  which  can  be  employed  to 
study  diffusive  as  We'Ll  as  hydrodynamic  flow  processes  when  occurring 
therein. 

The  model  proposed  in  this  work  satisfactorily  fulfils  the 
above  specified  objective;  it  permits  a  totally  rigorous  mathematical 
analysis  of  diffusive  as  well  as  hydrodynamic  flow  processes  and 
it  yields  encouraging  quantitative  predictions  in  both  instances. 
Moreover,  this  model  may  be  extended  to  the  study  of  transport 
processes  occurring  within  a  swarm  of  co-axially  orientated  spheroids. 
This  extension  has  permitted  an  assessment  of  the  important  effects 
which  particle  shape  and  orientation  can  have  in  systems  composed  of 
non-spherical  particles. 

The  solutions  to  some  important,  but  hitherto  unresolved, 
problems  have  incidentally  been  forthcoming  in  this  work,  of 
particular  interest  being  the  generalization  of  Stokes  Law  to  the 
case  of  the  porous  sphere.  In  addition,  an  exact  solution  to  the 
dual  problem  of  creeping  flow  through  a  porous  medium  containing  a 
spherical  cavity  has  been  effected. 
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NOMENCLATURE 


I.  SYMBOLS  PERTAINING  TO  DIFFUSIVE  FLOW  PROCESSES 


c 

D 

e 

K 

q 

Q 

R 

S 

[x,y,z] 
[r ,  0  ,<J>] 

[£,n,<j>] 


concentration  of  diffusing  species 

diffusivity  of  diffusing  species 

eccentricity  of  spheroid 

specific  conductivity 

diffusive  flux 

mainstream  diffusive  flux 

radius  of  sphere 

radius  of  unit  (model)  cell 

Cartesian  coordinates 

spherical  coordinates 

spheroidal  coordinates 


Greek  Symbols 

£  porosity  of  porous  medium 

t  tortuosity  of  pore  space 

X  diffusivity  (conductivity)  factor 

V  del  (gradient)  operator 

V2  Laplacian  operator 

Subscripts 

designates  vectorial  quantities 
designates  tensorial  quantities 
£  designates  randomly  orientated  spheroids 

s  designates  spheres  which  are  porous 

Superscript 

*  designates  macroscopically  averaged  quantities 

pertaining  specifically  to  a  porous  medium 
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II.  SYMBOLS  PERTAINING  TO  HYDRODYNAMIC  FLOW  PROCESSES 


F 

N 

P 

R 

S 

u 

U 

W 

[x,y,z] 

[r.e,*] 


hydrodynamic  resistance  force 
number  of  spheres  within  a  swarm 
pressure  referred  to  a  datum  plane 
radius  of  sphere 
radius  of  unit  (model)  cell 
velocity  of  flowing  fluid 
mainstream  velocity 

factor  by  which  the  resistance  of  a  swarm  of  spheres 
differs  from  that  predicted  by  Stokes  Law 
Cartesian  coordinates 
spherical  coordinates 
spherical  harmonic  operator 


Greek  Symbols 

a  normalized  radius  of  reference  sphere,  R//k 

S  normalized  radius  of  unit  cell,  S//< 

e  porosity  of  porous  medium 

k  permeability  of  porous  medium 

U  viscosity  of  fluid 

p  density  of  fluid 

ip  streamfunction 

t  tangential  shear  stress 

X  normalized  radial  coordinate,  r//i< 

v  normalized  radius  of  porous  sphere,  R/ /k 

V  del  (gradient)  operator 

V2  Laplacian  operator 

Subscript 

designates  vectorial  quantities 

Superscript 

*  designates  macroscopically  averaged  quantities  pertaining 

specifically  to  a  porous  medium 
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INTRODUCTION 


Broadly  speaking,  the  principal  transport  processes  occurring 
within  porous  media  may  be  classified  as  follows: 

Molecular  Transport 
Ionic  Transport 
Charge  Transport 
Thermal  Transport 

Momentum  Transport  }  Hydrodynamic  Flow  Processes 

The  principal  difference  between  these  two  major  categories  is  that 
whilst  diffusive  flow  processes  are  generally  described  by  the  Diffusion 
Equation  (in  the  steady  state  by  Laplace's  Equation),  hydrodynamic  flow 
processes  are  governed  by  the  Navier-Stokes  Equation. 

In  engineering  practice  these  fundamental  transport  processes 
often  occur  simultaneously.  For  example,  the  process  of  liquid  evapor¬ 
ation  within  a  porous  medium  would  encompass  molecular  transport,  thermal 
transport  and  momentum  transport.  However,  in  order  to  be  able  to  make 
meaningful  predictions  concerning  such  complicated  phenomena  a  fuller 
understanding  must  firstly  be  developed  of  the  fundamental  transport 
processes  when  occurring  individually  within  porous  media. 

An  appreciation  of  diffusive  and  hydrodynamic  flow  processes 
occurring  within  porous  media  (such  as  catalyst  beds,  sand  beds,  column 
packings,  etc.)  is  of  considerable  importance  in  many  fields  of 
engineering.  For  example,  in  design  work  a  prior  knowledge  of  the 
permeability  of  various  porous  media  is  a  great  asset.  Equally 
important  is  a  knowledge  of  the  effective  diffusivity  of  molecules  or 
ions  diffusing  through  a  porous  medium,  as  well  as  a  knowledge  of  the 
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specific  conductivity  of  the  medium  when  saturated  with  a  conducting 
fluid . 

It  is  well  known  that  the  accurate  measurement  of  absolute 
dif fusivities  is  excessively  difficult  on  account  of  the  inevitable 
presence  of  interfering  convection.  However,  measurements  of  the 
effective  diffusivity  within  a  porous  medium  (for  example,  a  randomly 
packed  bed  of  spheres)  are  usually  simpler  to  execute  since  convection 
effects  can  be  arbitrarily  reduced  by  a  judicious  choice  of  particle 
size.  A  knowledge  of  the  relationship  between  this  effective  diffusivity 
and  the  absolute  diffusivity  in  unobstructed  space  therefore  permits  an 
indirect  determination  of  this  latter  quantity. 

However,  transport  processes  occurring  within  porous  media 
rank  amongst  the  most  complex  phenomena  encountered  in  engineering 
practice.  This  is  because  the  internal  geometry  of  most  porous  media 
is  not  sufficiently  well  understood  to  afford  more  than  a  loose 
statistical  description.  Exact  analytical  studies  are,  therefore, 
usually  out  of  the  question  and  a  mathematical  model  must  be  developed 
to  predict,  and  to  yield  physical  insight  into,  such  phenomena.  That 
is,  certain  approximations  regarding  the  internal  geometry  must  be 
introduced  such  that  a  mathematical  solution  can  be  effected. 

Several  models  have  been  proposed  and  presented  in  the  litera¬ 
ture  to  study  specific  transport  processes  occurring  within  specific 
porous  media.  However,  there  does  not  appear  to  be  available  at  present 
any  one  model  which  can  be  used  to  predict  successfully  diffusive  as 
well  as  hydrodynamic  flow  processes  within  even  the  simplest  homogeneous 
and  isotropic  medium,  viz.  within  a  randomly  packed  bed  of  spheres. 
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This  is  disappointing  because  such  models  are  introduced  only  to 
simplify  the  complex  pore  structure  of  the  medium  in  question  and, 
ideally ,  their  overall  success  should  not  depend  on  which  particular 
transport  process  is  being  studied.  The  specific  system  of  an 
homogeneous  and  isotropic  swarm  of  spheres  is  of  particular  interest 
because  many  porous  media,  particularly  those  of  a  granular  nature, 
may  be  well  approximated  thereby  and  such  systems  can  be  readily 
prepared  in  the  laboratory  for  experimental  investigation;  moreover, 
such  a  system  can  be  defined  by  a  minimum  of  parameters,  viz.  the 
porosity  and  the  particle  size  distribution. 

A  considerable  number  of  models  aimed  at  predicting  transport 
processes  within  porous  media  have  been  based  on  reasonably  well 
developed  theories  for  these  processes  when  occurring  in  capillaries. 
However,  all  such  capillary  models  are  inherently  anisotropic  in 
constitution  and  without  exception  their  final  predictions  incorporate 
at  least  one  adjustable  parameter  which  must  be  determined  by  experiment. 
In  this  sense,  therefore,  approaches  based  on  the  capillary  model  can 
only  be  partially  predictive  and  can  often  be  misleading.  Moreover, 
such  models  are  not  productive  in  the  study  of  diffusive  processes. 

On  account  of  the  extreme  geometric  simplifications  inherent  in  the 
capillary  model  this  representation  requires  but  a  minimum  of  mathe¬ 
matical  sophistication  during  its  application.  In  contrast,  the  model 
proposed  later  in  this  work  will  contain  significantly  fewer  simpli¬ 
fications;  this  will  necessarily  lead  to  far  more  complicated  mathematics. 

The  principal  objective  of  this  analysis  is  to  develop  an  improved 
generalized  model  for  an  homogeneous  and  isotropic  swarm  of  spheres 
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which  can  be  used  to  predict  successfully  diffusive  as  well  as  hydro- 
dynamic  flow  processes  occurring  therein.  A  secondary  objective  will 
be  to  extend  this  model  to  the  study  of  anisotropic  systems  composed 
of  non-spherical  particles. 
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TEE  PROPOSED  MODEL  FOR  AN  HOMOGENEOUS  AND 
ISOTROPIC  SWARM  OF  SPHERES 

Figure  1  depicts  an  unbounded,  homogeneous  and  isotropic  swarm 
of  solid  spheres  possessing  a  porosity  e  and  an  arbitrary  size  distri¬ 
bution.  Choosing  any  reference  sphere  (of  radius  R)  within  the  swarm 
it  is  postulated  that  this  sphere,  together  with  its  associated  more 
space  (here  approximated  by  a  concentric  annulus  having  an  outer 
radius  S)  sees  the  remainder  of  the  system  as  an  homogeneous  and 
isotropic  porous  mass  of  porosity  e;  this  is  illustrated  in  Figure  2. 

In  other  words  the  portion  of  the  complicated  pore  space 
fairly  attributable  to  our  reference  sphere  has  been  replaced  by  a 
uniform  annulus,  the  outer  radius  of  which  must  be  specified  such  that 
the  porosity  of  the  unit  cell  (comprising  the  reference  sphere  and  its 
associated  annulus)  is  identical  to  that  prevailing  throughout  the 
original  system.  This  necessitates  that 

S/R  =  (1-e)"1 /3  0  <  e  <  1.0  .  (1) 

This  relation  ensures  that  the  uniform  porosity  e  is  not  locally 
disturbed  by  the  modelling  procedure.  An  inherent  advantage  of  the 
proposed  model  is  that  the  macroscopically  homogeneous  and  isotropic 
characteristics  of  the  original  system  remain  unchanged,  i.e.  the 
modelled  system  is  indistinguishable  from  the  original  system  when 
viewed  from  the  outside. 

To  summarize  then,  the  modelled  system  consists  of  our  reference 
sphere  and  two  adjoining  concentric  regions,  viz.  an  annular  region  of 
void  space  and  an  exterior  region  of  homogeneous  and  isotropic  porous 


ad;  .j.K  w  Hit  '  i  lania  10  **  1  >  a-i  a  ^ 


The  Proposed  Model 


6 


material.  In  essence,  the  technique  of  solution  will  be  to  solve  the 
differential  equations  which  describe  the  transport  process  of  interest 
within  each  of  these  regions  and  to  properly  connect  these  solutions 
at  the  interface  between  them.  It  is  important  to  note  that  subsequent 
to  the  above  indicated  modelling  of  the  porous  medium  the  mathematical 
developments  for  both  diffusion  and  fluid  flow  are  rigorous  in  their 
entirety. 
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FIGURE  1.  AN  HOMOGENEOUS  AND  ISOTROPIC  SWARM 
OF  SPHERES 

(Cross-Section  Through  Centre  of 
Reference  Sphere) 
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FIGURE  2.  THE  PROPOSED  MODEL  FOR  AN  HOMOGENEOUS 
AND  ISOTROPIC  SWARM  OF  SPHERES 
(Cross-Section  Through  Centre  of 
Reference  Sphere) 
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1.1  OBJECTIVES 

The  purpose  of  the  following  analysis  is  to  develop  a  fuller 
understanding  of  diffusion  taking  place  within  an  homogeneous  and  iso¬ 
tropic  porous  medium  composed  of  spherical  particles  possessing  an 
arbitrary  size  distribution;  in  particular  to  predict  the  diffusivity 
factor ,  X,  which  is  defined  to  be: 

Effective  diffusivity  within  porous  medium  D* 

X  =  - - - : - .  (2) 

Absolute  diffusivity  in  unobstructed  space  D 

This  will  be  accomplished  forthwith  by  application  of  the  proposed  cell 
model . 
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1.2  THE  DEFINING  EQUATIONS  AND 
BOUNDARY  CONDITIONS 

1.2.1  FUNDAMENTAL  DIFFERENTIAL  EQUATIONS 

The  modelled  system  for  diffusion  through  an  homogeneous  swarm 

of  spheres  is  depicted  in  Figure  3.  The  validity  of  Fick’s  Law  of 
2 

diffusion  within  the  annular  region  of  this  system  will  be  acknow¬ 
ledged,  viz. 

q=-DVc  R<r<S.  (3) 

2  27 

Furthermore,  Fick's  Law  in  its  macroscopic  form  ’  will  be  postulated 
as  a  description  of  conditions  within  the  exterior  region,  viz. 

q*  =  -D*  Vc*  S  <  r  <  00  .  (4) 


The  symbol  *  will  be  employed  throughout  this  study  to  designate 
macroscopically  averaged  quantities  pertaining  specifically  to  a 
porous  medium. 

In  the  steady  state  there  can  prevail  no  net  build-up  of  the 
diffusing  material;  this  implies  that 


V 


•  q  = 


0  , 


(5) 


V.q*  =  0  . 


(6) 


Since  D  and  D*  are  independent  of  location,  the  application  of  these 
continuity  conditions  to  Equations  (3)  and  (4)  yields: 

annular  region:  V2c  =0  R  <  r  <  S  ,  (7) 


exterior  region: 


V2c*  =  0 


S  <  r  <  00  . 


(8) 


. 
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FIGURE  3.  DESCRIPTION  OF  THE  MODELLED  SYSTEM  FOR  DIFFUSION  THROUGH  AN  HOMOGENEOUS 
SWARM  OF  SPHERES  (Cross-Section  Through  Centre  of  Reference  Sphere) 
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In  order  to  determine  the  global  concentration  distribution,  solutions 

of  Equations  (7)  and  (8)  are  sought  which  satisfy  physically  rational 

boundary  conditions  to  be  stipulated  later.  The  system  of  coordinates 

convenient  to  this  analysis  will  be  spherical  coordinates  £r,©,4>]  ,  in 

2  22 

which  the  Laplacian  operator  is  given  by  ’  : 

V2  =  (32/3r2)+(2/r)(3/3r)+(l/r2)(32/302)+(cot0/r2)(3/30) 

+(l/r2sin20)  (32/3(f>2)  ,  (9) 


and  Fick's  Laws  by: 


q  =  [q^q^q^]  =  [-D(3c/3r)  ,-D(l/r)(3c/36)  ,-D(l/rsin0)  (3c/3<J>)  ]  , 

(10) 

q*  =  [q*,q^,q|]  =  [-D*(3c*/3r)  ,-D*(l/r)(3c*/30)  ,-D*(l/rsin0)  (3c*/3<jO  ]  . 

(ID 

However,  for  macroscopically  rectilinear  diffusion  in  the  x-direction 
there  is  no  ^-dependence  so  that  =  0  and  q*  =  0.  Consequently  the 
<j)  variable  may  hereafter  be  suppressed. 

1.2.2  STIPULATED  BOUNDARY  CONDITIONS 

It  will  be  assumed  that  there  is  no  deposition  of  material  on 
the  reference  sphere  and  no  dissolution  thereof;  this  implies  that 

qrCR+,e>  =  o  .  (12) 

From  considerations  of  equilibrium  and  continuity  at  the  interface 
separating  the  annular  and  exterior  regions  it  is  necessary  that  the 
concentration  and  radial  flux  distributions  within  these  regions  match 


. 


- 
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identically  in  the  limit,  viz. 


c(s  ,e) 

=  c*(s+,e)  , 

(13) 

r(s\e) 

=  qJ(s+,0)  . 

(14) 

Moreover,  the  flux  vector  at  any  station  far  removed  from  the  reference 
sphere  must  approach  that  of  the  mainstream,  Q* ,  viz. 

Limit  q*(r,0)  =  Q*  ,  (15) 

r  ->  oo  — 

where  Q*  =  [Q*,0,0]  in  Cartesian  coordinates  [x,y,z]. 

It  is  not  possible  to  stipulate  any  boundary  conditions  concer¬ 
ning  q  and  q*  (the  tangential  components  of  the  flux  vectors)  owing 

0  U 

to  the  poorly  understood  phenomenon  of  surface  diffusion,  that  is  the 

9  27 

migration  of  molecules  or  ions  over  the  solid  surfaces  ’  .  It  is 

certainly  dangerous  to  write  q  (R,0)  =  0  by  direct  analogy  with  the 

w 

viscous  flow  problem  because  viscous  ftow  is  a  relative  motion  of  the 
adjacent  elements  of  a  fluid  whereas  diffusive  flow  is  a  relative 
motion  of  its  several  constituents. 

However,  it  transpires  that  Equations  (12)  -  (15)  constitute  an 
exact  set  of  physically  realistic  and  consistent  boundary  conditions, 
thereby  permitting  unique  solutions  of  Equations  (7)  and  (8) .  Subsequen¬ 
tly,  the  tangential  component  q  (R,0)  may  be  calculated  and,  moreover, 
it  does  purport  to  the  existence  of  surface  diffusion  over  the  solid 
surfaces  (see  Equation  (27)  to  follow) . 


« 


. 
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1.3  SOLUTION  OF  THE  DEFINING  EQUATIONS 

1.3.1  THE  GENERAL  SOLUTIONS 

The  general  solutions  of  Equations  (7)  and  (8)  for  the  case  of 

rectilinear  diffusion  may  be  extracted  from  fundamental  solutions  of 

22 

Laplace’s  Equation  in  spherical  coordinates  discussed  by  Sneddon  ,  viz. 

c(r,0)  =  £  [A  r  +  B  r  n  ■*■]  [P  (cos0)  +  C  Q  (cos0)  ]  R  <  r  <  S  ,(16) 
n 

c*(r ,0)  =  I[A*rn  +  B*r~n_1nP  (cos0)  +  C*  Q  (cos0)]  S  <  r  <  «  ,(17) 
n 

where  A  ,  B  ,  C  ,  A* ,  B*,  C*  represent  arbitrary  constants  and  P  (cos0) 
n  n  n  ri  it  n  n 

and  Q  (cos0)  denote  nth  order  Legendre  functions  of  the  1st  and  2nd 
n 

kind  respectively.  The  particular  solutions  which  satisfy  the  four 
stipulated  boundary  conditions  can  then  be  shown  to  be: 

c(r,6)  =  c*(r,7r/2)  +  A(RQ*/D*)  [  (r/R)+(l/2)  (r/R)  "2  ]  cos0  ,  (18) 

o*(r,0)  =  c*(r,ir/2)  +  (RQ*/D*)  [  (r/R)+(A*/2)  (r/R)_2]cos6  ,  (19) 

where  A  and  A*  are  given  by : 

A  =  3X/(3X-Xe+e)  ,  (20) 

A*  =  (3X-Xe-2e)/{(l-e)(3X-Xe+e)}t.  (21) 

In  developing  these  composite  expressions  for  A  and  A*  it  was  necessary 
to  substitute  for  the  ratio  (S/R)  according  to  Equation  (1). 

t  Refer  to  Section  1.3.4  for  further  information  regarding  the 
numerical  value  of  A*. 


'  -• 
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1.3.2  QUANTITATIVE  CONSISTENCY  OF  THE  MODELLED  SYSTEM 

The  solution  developed  above  represents  the  exact  solution  for 
the  physically  realizeable  system  of  a  single  sphere  concentrically 
suspended  inside  a  spherical  cavity  within  an  isotropic  porous  medium. 
Therefore,  if  this  simpler  system  is  to  be  quantitatively  represen¬ 
tative  of  the  original  swarm  of  spheres  it  is  necessary  that  within 
the  unit  cell  the  average  flux  in  the  mainstream  direction  be  equal 
to  that  of  the  mainstream  itself.  In  particular,  this  implies  that 

S 

(1/ttS2)  /  qQ(r,7T/2)  2nr  dr  =  Q*  .  (22) 

R 

The  flux  component  qQ(r,Tr/2)  may  be  evaluated  from  Equations  (10),  (18) 

O 

and  (20),  thus: 

qfi(r ,tt/2)  **  {-D(l/r)  (3c/30)  }  =  3Q*{l+(l/2)  (r/R) “3 }/(3A-Ae+e)  .(23) 

y  @tt/2 

Evaluating  the  integral  and  substituting  for  (S/R)  from  Equation  (1) 
yields  the  ultimate  expression 

A  =  D*/D  =  2e/(3-e)  0  <  e  <  1.0  .  (24) 

1.3.3  PHYSICAL  INSIGHT  PROVIDED  BY  THE  SOLUTION 

It  is  immediately  apparent  that  the  predictions  of  Equation  (24) 
are  physically  consistent  at  both  porosity  limits,  viz.  as  e  ->  0 , 

A  •*  0  and  as  e  ->  1,  A  ->  1.  Furthermore,  this  equation  implies  that  A 
is  invariant  with  the  size  distribution  for  systems  composed  entirely 
of  spheres  and  possessing  a  macroscopically  uniform  spatial  distribution. 
In  view  of  this  fact  it  is  clear  that  diffusivity  measurements  should 


. 


. 
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not  be  expected  to  provide  quantitative  information  relating  to  pore 
sizes,  pore  size  distributions  and  consequently  to  specific  surface 
areas  of  porous  media  in  general. 

1.3.4  DISTURBANCE  INTRODUCED  BY  TEE  MODELLING  PROCEDURE 

It  is  particularly  elucidating  to  calculate  the  extent  of  any 
disturbance  of  the  mainstream  flowfield  introduced  by  the  modelling 
procedure . 

Firstly,  from  Equations  (21)  and  (24)  it  may  in  fact  be  noted 

that  A *  =  0.  In  view  of  this  it  follows  from  Equations  (11)  and  (19) 

that  q*(r,0)  =  -Q^cosB  and  q*(r,0)  =  Q*sin0  .  The  flux  vector, 
r  0 

q*(r,0),  at  any  station  within  the  exterior  region  will  therefore  be 
characterized  by: 

|q*(r,0)|  =  /{ q* ( r , 0)  }2  +  (q*(r,0)}2  =  Q*  .  (25) 

Hence,  the  flowfield  everywhere  within  the  exterior  region  remains 
totally  undisturbed,  by  the  modelling  procedure.  In  other  words  any 
disturbance  introduced  by  the  modelling  procedure  is  wholly  confined 
to  within  the  unit  cell,  an  observation  which  lends  heavy  support  to 
the  realistic  nature  of  the  proposed  model. 

1.3.  5  EVALUATION  OF  TEE  SURFACE  DIFFUSION 

As  noted  previously  the  formulation  of  boundary  conditions 
concerning  the  tangential  components  of  the  flux  vector  was  avoided 
owing  to  the  insufficiently  understood  phenomenon  of  surface  diffusion. 
However,  the  magnitude  of  the  surface  flux,  q0(R,0),  over  the  reference 
sphere  may  now  be  evaluated  by  means  of  Equations  (10)  and  (18) ,  viz. 


. 


. 
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q0(R,6)  =  { — D( 1/ r) (3c/30) }  =  (9/2)Q*sin6/(3X-Xe+e)  .  (26) 

@R 

Substituting  for  X  from  Equation  (24)  then  yields  the  ultimate 
expression 

qe(R,6)/Q*  =  (3/ 2) sinG/ e  .  (27) 

This  result  suggests,  as  intuitively  expected,  that  the  contribution 
of  surface  diffusion  to  the  mainstream  flux  will  be  most  significant 
for  very  low  porosities. 

1.3.6  CONSIDERATIONS  WHEN  THE  SPHERES  THEMSELVES  ARE  POROUS 

Frequently  the  system  under  consideration  is  composed  of  porous 

particles,  for  example  a  bed  of  spherical  catalyst  pellets  in  which 

each  pellet  possesses  the  same  porosity,  e  ,  and  the  same  diffusivity 

s 

factor,  X  .  For  such  a  system  the  model  representation  assumes  the 
form  depicted  in  Figure  4.  To  ensure  that  the  porosity  of  the  unit 
cell  again  remains  equal  to  that  of  the  original  system,  e,  it  is  here 
necessary  that 

S/R  =  { (1-e) / (1-e  )}_1/3  .  (28) 

s 

The  calculations  for  this  system,  although  much  more  complicated, 
are  entirely  analogous  to  those  of  the  preceding  derivation  for  solid 
spheres.  However,  it  now  becomes  necessary  to  solve  Laplace's  equation 
within  three  regions,  viz.  the  exterior  region,  the  annular  region 
and  the  reference  sphere  itself,  and  to  connect  these  solutions  at  the 
appropriate  interfaces  by  means  of  boundary  conditions  analogous  to 
those  stipulated  in  Equations  (12)  -  (15).  The  ultimate  result  is: 


. 


. 


- 
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FIGURE  A. 


) 


-1/3 


THE  PROPOSED  MODEL  FOR  AN  HOMOGENEOUS  SWARM  OF 
POROUS  SPHERES 
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(2+A  )  -  2 ( 1—  X  )(l-e)/(l-e  ) 

A  =  - ® - i - i-  *  (29) 

(2+A  )  +  (1~A  ) (1-e) / (1-e  ) 

So  S 

For  the  particular  case  of  diffusion  through  solid  spheres 

(e  =  0,  A  =  0)  the  above  expression  reduces  to: 
s  s 

A  =  2e/(3-e)  .  (30) 

This  is  identical  with  Equation  (24)  for  solid  spheres  as  would  be 
expected . 

If,  moreover,  each  individual  sphere  may  be  visualized  as  being 

composed  of  considerably  smaller  spheres  then  A  may  be  evaluated 

s 

from  Equation  (30)  above,  viz. 

A  =  2e  / (3-e  )  .  (31) 

s  s  s 

The  substitution  of  this  expression  into  Equation  (29)  for  porous 
spheres  surprisingly  yields  the  original  expression  for  solid  spheres, 
viz . 

A  =  2e/(3-e)  .  (32) 

This  analysis  therefore  suggests  that  A  is  independent  of  e 

s 

for  the  special  case  of  porous  spheres  themselves  composed  of  much 
smaller  spheres.  Once  again  the  inference  is  that  A  is  invariant 
with  the  size  distribution  for  any  system  composed  entirely  of  spheres 
and  possessing  a  'macroscopically ?  uniform  spatial  distribution. 
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1.3.7  ELECTRICAL  CONDUCTION  THROUGH  AN  HOMOGENEOUS  SWARM  OF  SPHERES 

The  electrical  analogue  of  Equation  (32)  for  conduction  through 
an  homogeneous  swarm  of  solid,  non-conducting  spheres,  saturated  with 
a  conducting  fluid  of  specific  conductivity  K,  is 

A  =  K*/K  =  2e/(3-e)  ,  (33) 

where  A  here  denotes  the  conductivity  factor  of  the  porous  medium. 

That  this  identification  between  diffusion  and  electrical  conduction 
can  be  made  is  formally  apparent  from  the  similarity  in  structure  of 
the  underlying  boundary  value  problems,  namely  the  partial  differen¬ 
tial  equations  (based  upon  Fick's  Law  and  Ohm's  Law  respectively) 

and  the  stipulated  boundary  conditions^ .  Moreover,  Schofield  and 
.  19 

Dakshinamurti  have  actually  measured  the  diffusivity  and  conduct¬ 
ivity  factors  of  a  wide  range  of  sands  and  clays  and  have  obtained 
agreement  to  within  the  limits  of  experimental  error  (Klinkenberg^ 
has  further  confirmed  this  equivalence  from  similar  observations). 
Throughout  this  work,  therefore,  A  will  be  understood  to  be  synonomous 
with  the  diffusivity  factor  and  the  electrical  conductivity  factor. 


t  The  boundary  conditions  stipulated  for  diffusion  may  be 
interpreted  as  follows  for  electrical  conduction.  Thus,  Equation  (12) 
expresses  the  non-conducting  characteristics  of  the  reference  sphere, 
Equations  (13)  and  (14)  imply  respectively  potential  continuity  and 
radial-current  continuity  across  the  outer  surface  of  the  unit  cell, 
and  Equation  (15)  again  reflects  the  finite  nature  of  the  disturbance 
created  by  the  modelling  procedure. 

It  should  be  emphasized  that  the  expression  A  =  2e/(3-e)  will 
not,  in  general,  apply  to  the  case  of  thermal  conduction  because  no 
allowance  has  been  made  in  the  presented  derivation  for  either  a 
finite  film  resistance  between  the  spheres  and  the  fluid,  or  for  the 
inevitable  presence  of  convection. 


. 
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1.3.8  COMPARISON  OF  PREDICTED  RESULTS  WITH  EXPERIMENTAL  DATA 

There  exists  a  great  deal  of  electrical  conductivity  data  *  * 
supporting  Equation  (33)  for  e  >  0.4;  for  lower  porosities  the 
predictions  of  this  equation  are  somewhat  higher  than  observed  values 
(Figure  6) . 

As  regards  the  effects  of  the  size  distribution  of  the  component 

1  2_2.  26 

spheres,  the  consensus  of  opinion  *  *  is  that  X  is  essentially 

invariant  with  this  parameter  provided  the  spatial  distribution  of  the 

spheres  is  macroscopically  uniform;  these  observations  are  in  accord 

with  the  conclusions  of  the  present  analysis. 

However,  it  is  appropriate  at  this  juncture  to  stress  that 

the  experimental  measurement  of  the  conductivity  factor,  although 

significantly  easier  than  that  of  the  diffusivity  factor,  is  by  no 

means  trivial.  Such  electrical  measurements  are  often  hampered  by 

polarization  at  the  electrodes^  whilst,  for  dispersions  in  particular, 

certain  additional  circumstances  can  interfere  with  the  homogeneity 
14 

of  the  system  ,  notable  amongst  these  being  sedimentation  of  the 

dispersed  phase,  adherence  of  the  dispersed  phase  and  chain  formation 

in  the  presence  of  electric  fields.  The  very  small  effects  attributed 

14 

by  certain  workers  to  the  size  distribution  in  systems  of  spheres 
should,  therefore,  be  viewed  with  some  caution. 

Present  Experimental  Work 

In  order  to  confirm  the  trend  of  existing  experimental  data 
in  the  lower  porosity  range,  a  number  of  conductivity  factor  measure¬ 
ments  were  carried  out  using  a  cylindrical  plexiglass  cell  packed 


. 


■ 
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with  various  arbitrary  mixtures  of  glass  spheres  (Figure  5) .  Solid 

copper  electrodes  were  placed  at  both  ends  of  the  cell  and  an  acidified 

aqueous  copper  sulphate  solution  was  employed  as  the  electrolyte  in 

13 

order  to  effectively  eliminate  surface  effects  at  these  electrodes 

All  subsequent  electrical  measurements  were  made  using  audio-frequency 

alternating  current  in  order  to  avoid  polarization  effects.  In  order 

to  further  avoid  any  effects  due  to  packing  discrepancies  at  the  cell 

walls  it  was  necessary  to  ensure  a  minimum  cell  to  particle  diameter 

26 

of  at  least  twenty-five  to  one 

The  electrical  resistance  of  the  cell  when  filled  with  the 
electrolyte  alone  was  firstly  measured,  followed  by  the  resistance  of 
the  cell  when  packed  with  spheres  and  saturated  with  this  same  electro¬ 
lyte,  both  measurements  being  taken  at  the  same  temperature.  The  ratio 
of  these  two  resistances  gave  the  conductivity  factor  directly.  The 
porosity  of  the  system  was  finally  determined  by  a  weighing  technique. 
Representative  data  so  obtained  is  tabulated  below. 

TABLE  1 

EXPERIMENTAL  CONDUCTIVITY  FACTOR  DATA 
FOR  HOMOGENEOUS  SWARMS  OF  SPHERES 


e 

X 

0.261 

0.174 

0.296 

0.198 

0.309 

0.214 

0.348 

0.248 

0.380 

0.278 

This  data  is  also  displayed  in  Figure  6;  it  can  be  seen  to  conform 
closely  with  the  trend  of  previously  reported  experimental  data. 


' 
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C  CYLINDRICAL  PLEXIGLASS  CELL 
E  SOLID  COPPER  ELECTRODES 

S  PACK  OF  GLASS  SPHERES  SATURATED  WITH  ELECTROLYTE 
H  HOLE  FOR  LIQUID  EXPANSION 
R  STANDARD  RESISTANCE  100TI 

G  AUDIO  OSCILLATOR  1000Hz  (General  Radio  Co.  type  1311-A) 

D  DIGITAL  MULTI -FUNCTION  METER 
(Hewlett  Packard,  model  3450A) 


FIGURE  5.  EQUIPMENT  USED  FOR  THE  EXPERIMENTAL  DETERMINATION 
OF  THE  CONDUCTIVITY  FACTOR  OF  A  PACK  OF  SPHERES 
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FIGURE  6. 


THE  PREDICTED  DEPENDENCE  OF  A  ON  e  FOR  SPHERES: 
COMPARISON  WITH  EXPERIMENTAL  DATA 
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1.4  DISCUSSION  OF  PREVIOUS  WORK 

12 

Clerk  Maxwell  investigated  the  problem  of  electrical 
conduction  through  a  dispersion  of  solid  conducting  spheres  embedded 
within  a  conducting  medium.  In  essence  he  considered  the  limiting  case 
of  a  single  sphere  in  space  (e  1.0)  and  proposed  the  following 
solution  as  being  valid  for  dilute  systems  (e  -  1.0) ,  viz. 

(2+A  )  -  2 (1-A  ) (1-e) 

A  =  - - - - -  •  (34) 

(2+A  )  +  (1-A  ) (1-e) 

s  s 

This  equation  agrees  identically  with  the  electrical  conductivity 

analogue  of  the  derived  Equation  (29)  for  the  specific  case  of  solid 

conducting  spheres  (e  =  0,A  ^  0)  which  do  not  touch  one  another. 

s  s 

For  solid  non-conducting  spheres  (e  =  0,A  =0)  both  Equations  (29) 

s  s 

and  (34)  reduce  to  the  familiar  expression  derived  earlier,  viz. 

A  =  2e/(3-e)  .  (35) 

The  above  solution  is  closely  related  in  form  to  Lord  Rayleigh's"^ 

exact  solution  for  monosized  spheres  arranged  in  a  cubic  pattern,  viz. 

2e-0 . 3919 (1-e) 10 / 3  - 

A  =  - — - —  e  >  0.4764.(36) 

3-e-0. 3919 (1-e)  10/3  - - - - 

This  result  is  in  good  agreement  with  the  derived  Equation  (35) , 
differing  nowhere  by  more  than  3.02%  (Figure  7).  Lord  Rayleigh 
himself  questioned  whether  the  orientation  of  the  spheres  could  ap¬ 
preciably  affect  A  and  stated  that  "an  irregular  isotropic  arrangement 
would,  doubtless,  give  the  same  result". 
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Slawinski  studied  the  relationship  between  X  and  e  for 
aggregates  of  non-conducting  monosized  spheres  in  several  lattice 
orders.  However,  his  approach  was  based  on  merely  geometric  concepts 
and  thereby  disregarded  the  fundamental  principles  of  potential 
theory  according  to  which  the  potential  distribution  throughout  the 
conducting  medium  must  satisfy  Laplace's  Equation  (7).  Despite  this 
deficiency  Slawinski  was  able  to  develop  an  expression  which  is  in  good 
agreement  with  experimental  data  for  e  <  0.4,  viz. 

X  =  e/(l. 3219-0. 3219e)2  .  (37) 

3 

Bruggemann  examined  the  system  in  which  one  relatively  large 
sphere  is  surrounded  by  a  swarm  of  much  smaller  spheres.  He  considered 
the  region  exterior  to  this  large  sphere  to  be  a  continuum  and  subse¬ 
quently  applied  Maxwell's  result  on  the  premise  that  the  system  is 
'dilute'  with  respect  to  the  large  sphere.  He  derived  the  expression 

X  =  e3/2  .  (38) 


However,  for  single  sizes  or  narrow  size  fractions  the  physical 
conditions  necessary  for  justifying  the  Bruggemann  approximation  are 
not  satisfied,  and  this  is  reflected  in  the  less  satisfactory  over¬ 
all  agreement  of  his  formula  with  experimental  data  (Figure  7) . 

Archie^-  examined  contemporary  conductivity  data  for  spheres  and 
unconsolidated  sands  and  suggested  the  following  correlation: 

1  •  3 


X  =  £ 


(39) 


. 


■  l  .  >1  1  •  5 
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This  empirical  result  represents  the  data  well  for  e  <  0.4.  It  is 
interesting  to  note  the  similarity  in  form  of  Archie’s  Equation  (39) 
and  Bruggemann’s  Equation  (38). 

Of  particular  interest  is  the  work  of  Prager^  who  applied 
the  principle  of  minimum  entropy  generation  to  obtain  bounds  on  the 
diffusion  factor  for  an  homogeneous  suspension  of  solid  particles  of 
arbitrary  shape.  He  showed  that 

X  <  e{l-(l-e)/3}  ,  (40) 

and  stated  this  inequality  to  be  valid  for  particles  of  any  shape  and 
at  any  porosity,  the  only  stipulation  being  that  the  suspension  be 
isotropic.  For  spheres  in  particular  he  suggested  that 

X  =  e{l-(l-e)/2>  .  (41) 

This  formula  gives  good  agreement  with  experimental  data  over  the 
entire  porosity  range  (Figure  7) . 

It  is  instructive  to  note  that  Equation  (41)  constitutes  the 
first  two  terms  of  each  series  obtained  when  Equations  (38)  and  (35) 
are  expanded  in  terms  of  (1-e) ,  viz. 

X  =  £3/2  =  e{l-(l-e) l1 /2  =  e{l- (1-e) /2-(l-e) 2/8-(l-e) 3/16-* • • }  ,(42) 

X  -  2e/(3-e)  =  e{ l+(l-e) /2}~l  =  e{l-(l-e) /2+(l-e) 2/4-(l-e) 3/8+- • • }  . 

(43) 

These  expressions  therefore  exhibit  an  increasing  agreement  as  e  *>  1.0 
(Figure  7) . 
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yielded  expressions  of  the  form 


n  =  1  or  2  . 


(44) 


The  parameter  x  denotes  the  so-called  tortuosity  of  the  pore  space  *  ; 

this  quantity  must  be  measured  experimentally.  Clearly,  such  approaches 
are  impractical  since  the  determination  of  x  requires  a  prior  measure¬ 
ment  of  X  itself. 


It  has  recently  been  brought  to  my  attention  that  a  cell-type 


model,  fundamentally  related  to  that  proposed  in  this  work,  has  been 
employed  by  Hashin^  to  study  the  conductive,  magnetic  and  elastic 
properties  of  polycrystalline  aggregates  and  bi-metallic  composites. 
This  serves  to  illustrate  the  fundamental  significance  of  the  proposed 
model  in  representing  both  porous  and  non-porous  heterogeneous  media. 


‘ 


*. 
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FIGURE  7. 


THE  DEPENDENCE  OF  A  ON  e  FOR  SPHERES:  COMPARISON 
WITH  PREVIOUS  WORK 
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I . 5  SUMMARY 

The  presented  results  demonstrate  that  the  proposed  model 
offers  a  satisfactory  representation  of,  and  provides  valuable 
insight  into,  diffusion  occurring  within  an  homogeneous  swarm  of 
spheres.  For  solid  spheres  possessing  an  arbitrary  size  distri¬ 
bution  the  diffusivity  factor  may  be  evaluated  by  the  formula 
X  =  2e/(3-s).  If  the  spheres  are  non-conducting,  the  electrical 
conductivity  factor  also  may  be  evaluated  by  this  formula.  These 
predictions  agree  well  with  experimental  data  for  e  >  0.4;  for 
lower  porosities  the  predictions  are  somewhat  high.  The  overall 
agreement,  however,  is  satisfactory  and  lends  heavy  support  to 
the  realistic  nature  of  the  proposed  model  and  to  the  acceptability 
of  the  assumptions  implicit  therein. 

This  analysis  suggests  that  X  is  invariant  with  the  size 
distribution  of  the  spheres,  inferring  that  diffusivity  and  conduc¬ 
tivity  measurements  should  not  be  expected  to  yield  quantitative 
information  relating  to  pore  sizes,  pore  size  distributions  or 
specific  surface  areas  of  porous  media  in  general;  these  conclusions 
are  also  in  accord  with  experimental  observations. 

In  view  of  these  encouraging  results  attempts  will  now  be  made 
to  apply  the  proposed  model  to  the  study  of  diffusion  (conduction) 
occurring  within  certain  anisotropic  porous  media.  This  will  provide 
valuable  insight  into  the  important  effects  which  particle  shape  and 
orientation  have  on  the  diffusivity  (conductivity)  factor  of  uncon¬ 


solidated  porous  media. 


' 
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PART  IB 


EXTENSION  OF  THE  PROPOSED  MODEL 
TO  THE  STUDY  OF  DIFFUSION  WITHIN 
ANISOTROPIC  UNCONSOLIDATED  POROUS  MEDIA 
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1.6  OBJECTIVES 


The  principal  objective  of  this  analysis  is  to  assess  the 
effects  which  particle  shape  and  orientation  have  on  the  diffusivity 
factor  of  unconsolidated  porous  media. 

In  practice  many  such  porous  media  are  composed  of  preferen¬ 
tially  orientated  particles  far  removed  in  shape  from  the  spherical 
geometry,  for  example  a  filter  cake  consisting  of  flattish  or 
fibrous  particles.  An  excellent  representation  of  such  an  anisotropic 
system  is  often  provided  by  an  homogeneous  swarm  of  co-axially 
orientated  oblate  or  prolate  spheroids  respectively.  This  analysis 
seeks  to  achieve  the  above  mentioned  objectives  by  extending  the 
geometric  model  proposed  in  Part  IA  for  spheres  to  the  case  of 
co-axially  orientated  spheroids. 

However,  before  proceeding  further  with  the  problem  of  interest 
it  will  firstly  be  elucidating  to  examine  the  nature  of  the  differen¬ 
tial  equations  which  describe  diffusion  within  anisotropic  porous  media 


in  general. 


1.7  THE  EQUATIONS  GOVERNING  DIFFUSION 
THROUGH  ANISOTROPIC  POROUS  MEDIA 
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1.7.1  POROUS  MEDIA  IN  GENERAL 

It  is  widely  acknowledged  that  the  appropriate  form  of  Fick's 
Law  describing  steady  state  diffusion  within  anisotropic  porous  media 
(corresponding  to  Equation  (4)  for  isotropic  media)  is 


q*  =  -D*  Vc*  ,  (45) 

18 

where  D*  denotes  a  symmetric,  second-order  tensor  ;  in  a  Cartesian 
frame  of  reference  [x,y,z]  the  associated  matrix  would  be  represented 


by : 


D*  D*  D* 
xx  xy  xz 

D*  D*  D* 
xy  yy  yz 

D*  D*  D* 
xz  yz  zz 


(46) 


Corresponding  to  the  diffusivity  factor  for  isotropic  media,  A,  it 
here  becomes  necessary  to  define  a  tensor,  A,  thus: 


X 

X 

X 

XX 

xy 

xz 

X  =  (1/D) d*  = 

X 

X 

X 

—  — 

xy 

yy 

yz 

X 

X 

X 

xz 

yz 

zz 

(47) 


in  which 


X. .  =  (1/D)D*.  .  (48) 


The  diffusivity  D  again  refers  to  that  within  unobstructed  space. 

The  matrix  representation  in  Equation  (47)  will  be  diagonalized  by 

selecting  a  Cartesian  frame  of  reference  which  is  collinear  with 

18 

the  three  principal  directions  of  the  porous  medium  in  question  ,  thus 


, 


•  V 
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A  0  0 

x 

0  A  0 

y 

0  0  A 

z 


(49) 


1.7.2  POROUS  MEDIA  COMPOSED  OF  CO -AXIALLY  ORIENTATED  SPHEROIDS 


Having  resolved  matters  concerning  the  nature  of  the  tensor 
A  it  is  now  in  order  to  return  to  the  original  problem  of  interest, 
viz. the  study  of  diffusion  within  an  homogeneous  swarm  of  co-axially 
orientated  spheroids.  In  such  a  system  one  preferred  direction  of 
symmetry  will  be  observed;  thus,  if  the  x-axis  is  taken  to  be 
collinear  with  this  preferred  direction  (i.e.  collinear  with  the  axes 
of  revolution  of  the  individual  spheroids) ,  then  it  follows  from 
symmetry  that 


A  =  A  ,  (50) 

z  y 

so  that  the  diffusivity  factor  matrix  here  assumes  the  more  specific 
form 

(51) 

Hence,  in  order  to  extract  meaningful  information  concerning  diffusion 
through  such  an  anisotropic  system  it  will  suffice  to  investigate  the 
following  two  principal  cases,  viz. 


A  0 
x 


0 

0 


A 

y 

0  A 


Case  1:  Diffusion  collinear  with  the  axes  of  revolution  of  the  spheroids 
Case  2:  Diffusion  orthogonal  to  the  axes  of  revolution  of  the  spheroids. 


■ 
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The  spheroid  appears  to  be  the  most  appropriate  generalized 

geometry  to  study  because  (1)  it  approximates  a  vast  range  of 

physically  important  shapes  (ranging  between  flat  circular  discs, 

spheres  and  circular  cylinders) ,  and  (2)  it  lends  itself  to  a 

rigorous  mathematical  analysis.  Moreover,  the  oblate  spheroid  is 

now  regarded  to  be  a  somewhat  idealized  shape  in  many  mass  transfer 
20 

operations 


. 
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1.8  DIFFUSION  THROUGH  AN  HOMOGENEOUS  SWARM  OF 
CO-AXIALLY  ORIENTATED  OBLATE  SPHEROIDS 

1.8.1  THE  DERIVED  SOLUTIONS 

An  oblate  spheroid  is  generated  when  an  ellipse  is  rotated  about 
its  minor  axis;  such  a  geometry  may  be  visualized  as  a  sphere  which  has 
undergone  a  flattening  process.  Figure  8  depicts  an  homogeneous  swarm 
of  oblate  spheroids  of  identical  eccentricity,  e,  possessing  an  arbi¬ 
trary  size  distribution  and  orientated  such  that  the  axis  of  revolution 
of  each  particle  is  collinear  with  the  x-direction.  As  before,  A 

x 

denotes  the  diffusivity  factor  of  the  medium  collinear  with  this 
principal  direction  whilst  A^  denotes  the  diffusivity  factor  in  any 
orthogonal  direction  (i.e.  in  the  y-z  plane). 

The  modelled  system  for  co-axially  orientated  oblate  spheroids 
(Appendix  A,  Figure  27)  is  closely  related  to  that  for  spheres  (Figure 
3) ;  however,  the  eccentricity  of  the  outer  surface  of  the  unit  cell 
will  now,  in  general,  be  different  from  that  of  the  reference  particle. 
The  calculations  involved  are  analagous  to  those  for  spheres  but  are 
necessarily  much  more  tedious.  As  detailed  in  Appendix  A  the  following 
important  results  are  obtained  for  co-axially  orientated  oblate  spheroids 


=  { g(5Q)-g(51)  }/{gUQ)-f  (5^ }  , 

(52) 

=  {h($0)-h(51)}/{h(C0)-g(51)}  , 

(53) 

in  which  the  functions  f(0 ,  g(£)  and  h(£)  are  given  by: 


•  1 
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FIGURE  8. 


►  X 
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AN  HOMOGENEOUS  SWARM  OF  CO-AXIALLY  ORIENTATED 
OBLATE  SPHEROIDS  (Cross-Section) 
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f(£)  =  l/sinh$  -  arccot (sinh£)  ,  (54) 

g(5)  =  sinh5/cosh2£  “  arccot (sinh^)  ,  (55) 

h(5)  =  2f(5)  -  g(0  .  (56) 

The  parameters  and  ^  appearing  in  Equations  (52)  and  (53)  are 
defined  by : 

=  areatanh(e)  0  <  e  <  1.0  ,  (57) 

sinh3^^  +  sinh£^  -  (cosh25Qsinh^(_)) /(1-e)  =  0  .  (58) 

1.8.2  COMPUTED  RES ULTS 

Figure  9  records  the  A  =  A  (e)  and  A  =  A  (e)  relationships 

xx  y  y 

predicted  by  the  preceding  equations  for  several  representative  systems 

of  co-axially  orientated  oblate  spheroids. 

For  any  specific  eccentricity,  e,  it  may  be  observed  that 

A  <  A  ,  also  that  A  is  far  more  dependent  on  e  than  is  A  ;  this  is 
x  y  x  y 

especially  significant  for  e  <  0.5,  that  is  for  flattish  particles. 

For  the  particular  case  defined  by  e  =  1.0  the  oblate  spheroids 
are  spherical  and  the  predictions  offered  here  are  in  accord  with  those 
developed  specifically  for  spheres  (Section  1.3),  viz. 

A  =  A  =  2e/(3-e)  .  (59) 

x  y 

For  the  specific  case  defined  by  e  ->  0 ,  that  is  for  a  co¬ 
axially  orientated  swarm  of  very  thin  circular  discs,  the  present 


predictions  are: 


.  "•  ■ 
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0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1.0 


FIGURE  9.  THE  PREDICTED  DEPENDENCE  OF  X  AND  A  ON  e  FOR 

x  y 

CO-AXIALLY  ORIENTATED  OBLATE  SPHEROIDS 
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Ax  0  ,  (60) 

->  e  .  (61) 


This  representation  of  a  porous  medium  constitutes,  in  the  limit,  a 
plate-type  model. 
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1.9  DIFFUSION  THROUGH  AN  HOMOGENEOUS  SWARM  OF 
CO-AXIALLY  ORIENTATED  PROLATE  SPHEROIDS 

1.9.1  THE  DERIVED  SOLUTIONS 

A  prolate  spheroid  is  generated  when  an  ellipse  is  rotated  about 
its  major  axis;  such  a  geometry  may  be  visualized  as  a  sphere  which 
has  undergone  an  elongating  process.  Figure  10  depicts  an  homogeneous 
swarm  of  prolate  spheroids  of  identical  eccentricity,  e,  possessing  an 
arbitrary  size  distribution  and  orientated  such  that  the  axis  of 
revolution  of  each  particle  is  collinear  with  the  x-direction. 

As  detailed  in  Appendix  B  the  following  results  are  obtained 
for  co-axially  orientated  prolate  spheroids;  they  can  be  seen  to  be 
closely  related  in  form  to  those  for  oblate  spheroids  presented  in 
Equations  (52)  -  (58),  viz. 


=  {G(50)-G(^1)}/{G($0)-F(C1)}  , 

(62) 

=  {H(C0)-H(51)>/{H(50)-G(C1)}  , 

(63) 

in  which  the  functions  F(0,  G(£)  and  H(£)  are  given  by: 


F(C)  =  l/cosh£  -  areacoth(coshC)  , 

G(£)  =  cosh£/sinh2£  -  areacoth(cosh£)  , 

H (O  =  2F(5)  -  G(0  . 

The  parameters  £q  and  ^  appearing  in  Equations  (62)  and  (63) 
defined  by: 

=  areatanh(e)  0  <  e  <  1.0  , 

cosh3^  -  coshS^  -  (sinh2C0cosh50) / (1-e)  =  0  . 


(64) 

(65) 

(66) 

are  here 


(67) 


(68) 


'  :  • 
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-►  X 


FIGURE  10.  AN  HOMOGENEOUS  SWARM  OF  CO-AXIALLY  ORIENTATED 
PROLATE  SPHEROIDS  (Cross-'Section) 
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1.9.2  COMPUTED  RESULTS 

Figure  11  records  the  A  =  A  (e)  and  A  =  A  (e)  relationships 

xx  y  y 

predicted  by  the  preceding  equations  for  several  representative 

systems  of  co-axially  orientated  prolate  spheroids.  For  any  specific 

eccentricity,  e,  it  may  be  observed  that  A  >  A  ;  however,  in  contrast 

x  y 

to  the  preceding  case  for  oblate  spheroids,  A  is  only  slightly  more 

X 

dependent  on  e  than  is  A  . 

y 

The  predictions  for  e  =  1.0,  that  is  for  spheres,  are  once 

again : 


A  =  A  =  2e/(3-e)  .  (69) 

x  y 

However,  for  e  ->  0 ,  that  is  for  a  co-axially  orientated  swarm 
of  very  thin  circular  cylinders  the  predictions  are  as  follows: 

A  ->  e  ,  (70) 

X 

Ay  ->  e/(2-e)  .  (71) 


This  representation  of  a  porous  medium  constitutes,  in  the  limit,  a 
capillary- type  model. 

Equation  (70)  does  of  course  represent  the  exact  solution  for 
diffusion  parallel  to  co-axially  orientated,  very  long  cylinders. 
Moreover,  Equation  (71)  constitutes  a  truncation  of  Lord  Rayleigh's^ 
exact  solution  for  conduction  perpendicular  to  monosized  circular 
cylinders  arranged  in  a  square  pattern,  viz. 


£  -0.3058(l-s)4-0.0134(l-£) 8 - 

y  2-£-0.3058(l-e)4-0.0134(l-e) 8 - 


£  >  0.2146  .  (72) 


' 


'  • 


. 
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FIGURE  11.  THE  PREDICTED  DEPENDENCE  OF  X  AND  X  ON  e  FOR 

x  y 

CO-AXIALLY  ORIENTATED  PROLATE  SPHEROIDS 
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This  prediction  is  displayed  in  Figure  11.  It  should  be  emphasized 
that  Rayleigh’s  solution  applies  only  to  monosized  cylinders  occupying 
a  square  array,  for  which  0.2146  represents  the  minimum  obtainable 
porosity.  However,  for  co-axially  orientated  cylinders  possessing  a 
wide  distribution  of  radii,  arbitrarily  lower  porosities  are  attainable. 


Practical  Significance  of  the  Predicted  Results 

It  is  appropriate  to  mention  here  that  knowledge  of  the  ratio 

A x/ provides  valuable  information  concerning  the  extent  of  anisotropy 

in  systems  composed  of  co-axially  orientated  oblate  or  prolate  spheroids. 

The  predicted  results  have  important  significance  because  it  sometimes 

becomes  necessary  to  actually  construct  a  porous  medium  possessing  a 

specified  degree  of  anisotropy  (for  example  a  specified  value  of 

A  /A  ).  Such  a  system  may  often  be  approximated  by  preparing  an 
x  y 

homogeneous  swarm  of  co-axially  orientated  disc-like  or  rod-like 
particles  (approximating  oblate  and  prolate  spheroids  respectively) , 
an  estimate  of  the  required  ’eccentricity'  of  these  particles  being 
obtained  by  inspection  of  Figure  9  (for  discs)  and  Figure  11  (for 


rods) . 


. 
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I. 10  DIFFUSION  THROUGH  A  SWARM  OF  CO-AXIALLY  ORIENTATED  SPHEROIDS 

AT  AN  ARBITRARY  ANGLE  OF  ATTACK 

The  law  describing  diffusion  through  an  anisotropic  porous 
medium  composed  of  co-axially  orientated  spheroids  was  discussed 
earlier  in  context  with  Equations  (45)  -  (51) .  In  the  preferred 
frame  of  reference  [x,y,z]  it  assumes  the  specific  form: 

Vc*  .  (73) 

It  is  clear  that  the  predictions  for  X  and  X  offered  in 

x  y 

Sections  1.8  and  1.9  may  be  used  to  construct  the  diffusivity  factor 
matrix  appearing  in  Equation  (73)  above.  Moreover,  this  matrix 
provides  a  knowledge  of  the  tensor  X  which,  in  turn,  permits  a  full 
description  of  the  diffusion  (conduction)  occurring  as  a  result  of 
a  concentration  (potential)  gradient  in  any  arbitrary  direction, 
within  any  chosen  frame  of  reference. 


q*  =  -D 


X  0 
x 


0 

0 


0  X 

y 

o  o  x 
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I. 11  DIFFUSION  THROUGH  AN  HOMOGENEOUS  SWARM 
OF  RANDOMLY  ORIENTATED  SPHEROIDS 

1. 11.1  DISCUSSION  OF  THE  SYSTEM 

Under  certain  circumstances  the  system  under  consideration  may 
be  better  approximated  by  a  swarm  of  randomly  orientated  spheroids 
(which  in  consequence  will  be  isotropic)  than  by  a  swarm  of  co-axially 
orientated  spheroids.  This  is  frequently  true  of  particulate  suspen¬ 
sions.  Seldom  would  it  apply  to  sedimented  media,  as  in  such  systems 
the  individual  particles  would  usually  have  settled  with  some  definite 
bias  as  regards  preferential  orientation. 

The  solution  for  such  a  system  of  randomly  orientated  spheroids 
may  be  constructed  by  a  weighted  superposition  of  the  three  principal 

solutions  for  the  co-axial  case.  Now,  by  comparing  Fick’s  Law  with 
2  -1 

Ohm’s  Law  it  becomes  apparent  that  X  is  a  direct  measure  of  the 
resistance  of  a  porous  medium  to  diffusion.  The  following  estimate 
may  therefore  be  made  for  the  diffusivity  (conductivity)  factor, 

X^, ,  of  a  system  of  randomly  orientated  spheroids  of  identical  eccen¬ 
tricity,  possessing  an  arbitrary  size  distribution: 

X"1  =  (1/3) X_1  +  (2/3) X_1  ,  (74) 

E  x  y 

where  X  and  X  are  defined  by  Equations  (52)  and  (53)  for  oblate 

x  y 

spheroids,  and  by  Equations  (62)  and  (63)  for  prolate  spheroids. 

1. 11 . 2  COMPUTED  RESULTS 

Figures  12  and  13  display  the  predicted  X  =  X  (e)  relationships 
for  several  representative  systems  of  randomly  orientated  oblate  and 
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FIGURE  12. 


THE  PREDICTED  DEPENDENCE  OF  X  ON  e  FOR  RANDOMLY 
ORIENTATED  OBLATE  SPHEROIDS 
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FIGURE  13. 


THE  PREDICTED  DEPENDENCE  OF  ^  ON  e  FOR  RANDOMLY 
ORIENTATED  PROLATE  SPHEROIDS 
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prolate  spheroids  respectively.  It  is  immediately  apparent  that  for 
oblate  spheroids  A  is  strongly  dependent  on  the  eccentricity,  whereas 

Li 

for  prolate  spheroids  this  dependence  is  surprisingly  mild. 

The  Limiting  Solutions  for  Oblate  Spheroids 
The  solution  for  e  =  1.0  (spheres)  follows  directly  from 
Equations  (59)  and  (74),  viz. 

A  =  2e/(3-e)  .  (75) 

Lx 

This  result  is  to  be  anticipated  because  spheres  do  not  possess  a 
preferred  axis  of  orientation. 

The  solution  for  e  ->  0  (very  thin  circular  discs)  follows 
from  Equations  (60),  (61)  and  (74),  viz. 

A  0  .  (76) 

This  result  is  to  be  expected  from  physical  considerations. 

The  Limiting  Solutions  for  Prolate  Spheroids 

The  solution  for  e  =  1.0  (spheres)  follows  from  Equations 
(69)  and  (74)  and  is  identical  to  the  prediction  of  Equation  (75) 
above,  viz. 

A  =  2c/ (3-e)  .  (77) 

Lx 

Again,  this  result  is  to  be  anticipated. 

The  solution  for  e  ->  0  (very  thin  circular  cylinders)  follows 
from  Equations  (70),  (71)  and  (74),  but  cannot  be  anticipated,  viz. 

A  ->  3c/  (5-2e)  .  (78) 

L» 

Xt  is  interesting  to  observe  the  similarity  in  structure  of 


Equations  (77)  and  (78) . 
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1.12  DISCUSSION  OF  PREVIOUS  WORK 

No  work  has  been  reported  to  date  concerning  diffusion  or 

conduction  through  a  swarm  of  co-axially  orientated  spheroids. 

5  12 

However,  Fricke  did  extend  Maxwell's  theory  (concerning  conduc¬ 
tion  through  a  dispersion  of  spheres)  to  the  case  of  randomly 
orientated  spheroids.  From  considerations  of  the  potential  of 
a  single  spheroid  in  space  he  developed  the  following  formula 
for  a  low  concentration  dispersion  of  non-conducting  spheroids 
within  a  conducting  medium: 

=  Ae/(A+l-e)  ,  (79) 

the  parameter  A  being  a  function  of  eccentricity  alone.  Some  specific 
values  of  A  for  different  shapes,  as  computed  from  Fricke 's  compli¬ 
cated  formulae,  are  tabulated  below  for  reference. 


TABLE  2 


PREDICTIONS  FOR  RANDOMLY  ORIENTATED  SPHEROIDS 

(after  Fricke^) 


e 

A 

Ob  late 

Prolate 

1.0 

2.000 

2.000 

0.9 

1.994 

1.995 

0.8 

1.974 

1.979 

0.7 

1.931 

1.950 

0 . 6 

1.854 

1.909 

0.5 

1.730 

1.854 

0.4 

1.541 

1.786 

0.3 

1.271 

1.707 

0.2 

0.911 

1.623 

0.1 

0.474 

1.545 

0.0 

0.000 

1.500 

< 


. 


-Jr\  0  «i>LllHV 


, 


' 


1.12] 


53 


Fricke  himself  was  of  the  opinion  that  his  predictions  were 
rather  high,  considerably  so  for  oblate  spheroids^’7. 

It  is  interesting  to  note  that  the  predictions  offered  here 
agree  identically  with  those  of  Fricke  in  the  limit  e  -»■  1.0;  this 
agreement  is  to  be  anticipated  since  both  theories  then  effectively 
consider  a  single  spheroid  in  infinite  space.  The  present  predictions 
are  also  identical  with  Fricke' s  for  the  specific  geometries  defined 
by  e  =  1.0  (spheres)  and  e  •*  0  (thin  circular  discs  or  cylinders). 

For  all  other  values  of  eccentricity  the  present  predictions  (for  all 
porosities  0  <  e  <  1)  are  smaller  than  those  of  Fricke,  considerably 
smaller  for  oblate  spheroids. 

These  characteristics  of  the  presented  theory  appear  to  be 
a  direct  consequence  of  the  fact  that  the  eccentricity  of  the  outer 
surface  of  our  model  cell  differs  from  the  eccentricity  of  our  reference 
spheroid  (Appendix  A,  Figure  27)  for  all  geometries  excepting  those 
defined  by  e  =  1.0  and  e  ->  0 ;  for  these  specific  geometries  (viz.  spheres, 
thin  circular  discs  and  cylinders)  the  eccentricity  of  the  outer  sur¬ 
face  of  our  model  cell  is  identical  with  that  of  our  reference 
spheroid,  observations  which  follow  directly  from  Equations  (58)  and 
(68). 

Figure  14  displays  the  experimental  conductivity  data  of  De  La 
Rue  and  Tobias^  for  a  dispersion  of  randomly  orientated  cylinders 
possessing  a  diameter  to  length  ratio  of  0.1.  Now,  such  a  cylinder  is 
well  approximated  by  a  prolate  spheroid  possessing  an  eccentricity 
of  0.1;  the  results  predicted  by  Equation  (74)  for  randomly  orientated 
prolate  spheroids  possessing  this  same  eccentricity  can  be  seen  to 
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conform  closely  with  the  presented  data. 

It  will  be  recalled  that  the  predictions  of  Equation  (74) 
for  randomly  orientated  prolate  spheroids  lie  within  a  very  narrow 
band  over  the  entire  porosity  range  (Figure  13) .  That  the  data  of 
De  La  Rue  and  Tobias  for  highly  ’eccentric'  cylinders  lies  within 
this  band  lends  heavy  support  to  the  validity  of  the  proposed  theory. 
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FIGURE  14.  THE  PREDICTED  RESULTS  FOR  RANDOMLY  ORIENTATED 

PROLATE  SPHEROIDS:  COMPARISON  WITH  EXPERIMENTAL 
DATA 
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1.13  SUMMARY 

The  presented  results  demonstrate  that  the  proposed  cell 
model  offers  a  satisfactory  representation  not  only  of  an  isotropic 
swarm  of  spheres  but  also  of  an  important  class  of  anisotropic 
systems  (viz.  those  composed  of  co-axially  orientated  disc-like  or 
rod-like  particles) . 

In  particular,  it  has  been  demonstrated  that  swarms  of  co-axially 
orientated  spheroids  possessing  an  eccentricity  of  less  than  0.5  are 
markedly  anisotropic,  although  this  is  much  more  in  evidence  with 
oblate  than  with  prolate  spheroids.  These  results  have  important 
practical  significance  in  those  areas  in  which  conductivity  measure¬ 
ments  are  made  on  systems  composed  of  flatfish  or  fibrous  particles, 
a  notable  example  being  on  the  deliberately  deposited  mudcakes  in 
oil-well  bores. 

A  study  has  also  been  made  of  diffusion  (conduction)  occurring 
within  an  isotropic  system  composed  of  randomly  orientated  spheroids. 

The  diffusivity  (conductivity)  factor  here  has  been  demonstrated  to 
be  strongly  dependent  on  the  eccentricity  for  oblate  spheroids,  but 
surprisingly  invariant  with  this  parameter  for  prolate  spheroids. 

In  view  of  these  further  encouraging  results  attempts  will  now 
be  made  to  apply  the  proposed  cell  model  to  the  study  of  a  radically 
different,  and  far  more  complicated,  transport  process.  This  will 
involve  incompressible  creeping  flow  through  an  homogeneous  swarm 
of  spheres  and  will  provide  valuable  insight  into  the  nature  of 
liquid  flow  through  porous  media  in  general. 


- 
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PART  II 


STEADY  STATE  INCOMPRESSIBLE  CREEPING  FLOW 
THROUGH  AN  HOMOGENEOUS 
AND  ISOTROPIC  SWARM  OF  SPHERES 


, 
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II. 1  OBJECTIVES 

The  purpose  of  this  investigation  is  to  develop  a  fuller 
understanding  of  fluid  flow  occurring  within  an  homogeneous  and 
isotropic  porous  medium  composed  of  spherical  particles  possessing 
an  arbitrary  size  distribution;  in  particular,  to  evaluate  the 
resistance  offered  by  such  a  system  to  an  incompressible  fluid!  in 
the  creeping  flow  regime.  In  theory,  the  total  resistance  offered 
by  such  a  system  may  be  calculated  by  integrating  the  local  shear 
stress  over  the  entire  surface  thereof.  However,  so  excessively 
complicated  is  the  internal  geometry  of  even  the  most  regular  array 
of  spheres  that  resort  must  be  made  to  a  simplifying  model. 

On  account  of  the  encouraging  results  obtained  during  the 

preceding  diffusion  analyses,  the  proposed  cell  model  will  here 

be  further  employed  to  evaluate  the  permeability  of  an  homogeneous 

swarm  of  spheres  in  the  creeping  flow  regime.  For  liquid  flow  (with 

which  this  investigation  will  be  exclusively  concerned)  through 

porous  media  it  is  generally  acceptable  to  make  this  creeping  flow 
2  8 

assumption  ,  namely  that  inertial  forces  may  be  neglected  in 
comparison  with  viscous  forces. 

However,  before  proceeding  further  with  the  problem  of 
interest  it  will  firstly  be  elucidating  to  examine  the  nature  of  the 
equations  which  are  normally  employed  to  describe  liquid  flow  through 
porous  media,  and  which  implicitly  define  the  permeability  thereof. 

t  Throughout  this  work  attention  will  be  confined  to  Newtonian 
fluids  under  isothermal  conditions. 


- 
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II. 2  THE  EQUATIONS  GOVERNING  INCOMPRESSIBLE 
CREEPING  FLOW  THROUGH  POROUS  MEDIA 

There  has  been  considerable  dissension  in  the  past  concerning 
the  nature  of  the  differential  equation  which  describes  fluid  flow 
through  porous  media.  As  a  basis  for  further  discussion  it  will 
here  be  convenient  to  consider  the  somewhat  idealized  system  depicted 
in  Figure  15,  viz.  an  extent  of  homogeneous  and  isotropic  porous 
material  adjacent  to  a  region  of  unobstructed  fluid  space. 

Free  Fluid  Region 

Throughout  this  region  the  validity  of  the  Navier-S tokes 

Equation  will  be  acknowledged.  Thus,  for  the  flow  of  a  Newtonian 

fluid  at  a  sufficiently  low  Reynolds  number  the  hydrodynamic  con- 

28 

ditions  within  this  region  will  be  described  by  the  equation  : 

yV2u  =  Vp  .  (80) 

In  this  equation  u  denotes  the  velocity  vector,  p  the  absolute 
viscosity  of  the  fluid  and  p  the  pressure  referred  to  a  datum  plane. 

Porous  Medium  Region 

42 

The  following  hypothetical  differential  form  ,  based  on  the 
pioneering  experiments  of  Darcy  in  1851,  is  generally  quoted  with 
confidence  to  describe  incompressible  creeping  flow  within  an 
isotropic  porous  medium,  viz. 


-(p/<)_u*  =  Vp*  . 


(81) 
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FIGURE  15.  AN  IDEALIZED  TWO-REGION  SYSTEM 
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In  this  equation  k  denotes  the  pex’meab't'l'ity  of  the  porous  medium, 
u~  the  macroscopic  velocity  vector  and  p*  the  local  mean  pressure 
referred  to  the  same  datum  plane  (the  symbol  *  will  continue  to 
designate  macroscopically  averaged  quantities  pertaining  specifically 
to  a  porous  medium) .  The  velocity  profile  within  the  porous  medium 
region  according  to  this  equation  (now  universally  termed  the  Darcy 
Equation)  is  depicted  in  Figure  15. 

30 

As  long  ago  as  1949  Brinkman  proposed  that  a  viscous  stress 
term  need  be  incorporated  in  the  Darcy  Equation  (81)  in  order  to 
account  for  the  distortion  of  macroscopic  velocity  profiles  in  the 
vicinity  of  containing  walls  and  regions  of  unobstructed  fluid  space, 
thus  : 

-(p/k)u*  +  yV2u*  =  Vp*  .  (82) 

This  equation  (hereafter  referred  to  as  the  Brinkman  Equation)  has 

recently  received  encouraging  theoretical  substantiation  from 

54  55 

Slattery  and,  quite  independently,  from  Tam  .  For  small  values  of 

k  the  Brinkman  Equation  (82)  is  approximated  by  the  Darcy  Equation  (81) 
whilst  for  k  ->  00  it  possesses  the  inherent  advantage  of  developing  into 
the  Navier-S tokes  Equation  (80) ;  this  behaviour  is  to  be  anticipated 
because  a  non-porous  mass  and  an  unobstructed  fluid  represent  the 
possible  extremes  of  the  porous  medium.  The  velocity  profile  accord¬ 
ing  to  the  Brinkman  Equation  (82)  is  displayed  schematically  in 
Figure  15.  It  should  be  noted  that  this  equation  further  reduces  to 
the  Darcy  Equation  (81)  whenever  V2ju*  =  0,  that  is  whenever  the 
velocity  vector  _u*  is  uniformly  constant  throughout  the  porous  medium 
in  question  (as  has  been  presumed  in  Darcy’s  pioneering  experiments). 


- 
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Considerations  at  the  Interface 

At  the  interface  separating  the  two  regions  it  is  necessary 
that  the  limiting  solution  of  the  micro- flowfield  associated  with 
the  free  fluid  region  (as  described  by  the  Navier-Stokes  Equation) 
must  match  correctly,  when  considered  as  a  macro-f lowfield ,  with 
the  macro- flowfield  associated  with  the  porous  medium  region  (as 
described  by  either  the  Darcy  Equation  or  the  Brinkman  Equation) , 
thereby  reflecting  physical  consistency  (Figure  15) .  This  consti¬ 
tutes  the  principal  reason  as  to  why  the  Darcy  Equation  must  be 
regarded  as  being  incomplete  and  consequently  inadequate,  the 
explanation  being  that  the  Navier-Stokes  Equation  (80)  is  a 
differential  equation  of  the  second  order  whilst  the  Darcy  Equation 

(81)  is  of  the  first  order  only;  this  makes  it  impossible  to  formulate 
physically  rational  boundary  conditions  at  the  interface  using  the 
Darcy  representation  within  the  porous  region.  However,  this 
difficulty  is  not  encountered  when  employing  the  Brinkman  Equation 

(82)  ,  which  is  indeed  of  the  second  order  by  virtue  of  the  presence 
of  its  viscous  stress  term;  the  inclusion  of  such  a  term  is  therefore 
imperative  if  physical  consistency  is  to  be  maintained  throughout  the 
two  region  system  depicted  in  Figure  15 t . 


t  In  order  to  demonstrate  the  consistency  and  effectiveness  of 
the  Brinkman  Equation  (82)  three  problems  of  considerable  practical 
importance,  which  cannot  be  solved  rigorously  using  the  Darcy  Equation 
(81),  have  been  examined  in  some  detail  in  Section  II. 9,  Appendix  D 
and  Appendix  E.  These  problems  relate  respectively  to  (1)  sedimen¬ 
tation  of  an  isolated  porous  sphere,  (2)  incompressible  creeping  flow 
parallel  to  a  fracture  within  a  petroleum  reservoir,  and  (3)  incompres¬ 
sible  creeping  flow  through  an  isotropic  porous  medium  containing  a 
spherical  cavity. 
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Having  resolved  matters  concerning  the  nature  of  the  differential 
equation  which  governs  creeping  flow  through  porous  media  it  is  now  in 
order  to  return  to  the  original  problem  of  interest,  viz.  the  applica¬ 
tion  of  the  proposed  cell  model  (Figures  1  and  2)  to  the  study  of 
incompressible  creeping  flow  through  an  homogeneous  swarm  of  spheres. 
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II. 3  THE  MODELLED  SYSTEM  FOR  FLUID  FLOW  THROUGH 
AN  HOMOGENEOUS  SWARM  OF  SPHERES 

The  model  representation  for  fluid  flow  through  an  homogeneous 
swarm  of  spheres  is  depicted  in  Figure  16;  this  may  be  compared  with 
Figure  3  for  the  corresponding  problem  of  diffusion.  As  before,  it 
consists  of  a  reference  sphere ,  an  annular  region  of  void  space  and 
an  exterior  region  of  homogeneous  and  isotropic  porous  material. 

The  radius  of  the  unit  cell  (comprising  the  reference  sphere  and 
the  annular  region)  must  again  be  related  according  to  Equation  (1) , 
viz. 


S/R  =  (l-e)"1/3  0  <  e  <  1.0  ,  (83) 

which  ensures  that  the  uniform  porosity  e  is  not  locally  disturbed 
by  the  modelling  procedure. 
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FIGURE  16.  DESCRIPTION  OF  THE  MODELLED  SYSTEM  FOR  FLUID  FLOW  THROUGH  AN  HOMOGENEOUS 
SWARM  OF  SPHERES  (Cross-Section  Through  Centre  of  Reference  Sphere) 
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II. 4  THE  DEFINING  EQUATIONS  AND  BOUNDARY  CONDITIONS 

It  will  be  noted  during  the  development  which  follows  that 
the  fundamental  differential  equations  and  the  stipulated  boundary 
conditions  are  radically  different  from,  and  far  more  complicated  than, 
those  for  the  corresponding  diffusion  problem.  It  should  also  be 
stressed  that  the  subsequent  analysis  is  rigorous  in  its  entirety, 
requiring  no  physical  or  mathematical  simplifications. 

II. 4.1  FUNDAMENTAL  DIFFERENTIAL  EQUATIONS 

Within  the  annular  region  of  the  modelled  system  (Figure  16) 
the  Navier-S tokes  Equation  (80)  will  be  acknowledged  to  describe  the 
prevailing  creeping  flowfield,  thus: 

yV2u  =  Vp  R  <  r  <  S  .  (84) 

Within  the  exterior  region  the  validity  of  the  Brinkman  Equation 
(82)  will  be  proposed,  thus: 

-  (  y  /  k  )  u*  +  yV2u*  =  Vp*  S  <  r  <  «  .  (85) 

.  28 

For  steady  state  incompressible  flow  the  continuity  equation 
assumes  the  following  specific  forms  for  the  annular  and  exterior 
regions  respectively: 


V. u  =  0  , 

(86) 

V.  u*  =  0  . 

(87) 

The  terms  involving  pressure  in  Equations  (84)  and  (85)  may 
be  annihilated  by  performing  the  curl  vector  operation  throughout. 
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By  noting  that  pressure  is  a  scalar  quantity,  and  that  y  and  k  are 
independent  of  location,  the  following  expressions  are  obtained,  viz. 

y[VxV2u]  =  [VxVp]  =  0  ,  (88) 

~(u/k)[Vxu*]  +  p[VxV2u*]  =  [V_xVp*]  =  0  .  (89) 

The  system  of  Euler  (stationary)  coordinates  convenient  to 
this  investigation  will  be  spherical  coordinates  [r,0,cj)],  for  which: 


^ =  [VVud  - 

u*  =  [u*,u*,u*]  . 


(90) 

(91) 


The  expansions  of  the  individual  terms  in  Equations  (88)  and  (89)  in 

spherical  coordinates  are  presented  in  numerous  texts,  notably  Bird 
2  8 

et  al  .  On  account  of  the  specifications  introduced  so  far  the 

prevailing  flowfield  will  exhibit  axial  symmetry.  Thus,  u^  =  0  and 

u*  =  0;  hence,  the  <j>  variable  may  hereafter  be  supressed.  Consequently, 

streamf unctions  \p  and  ip*  may  now  be  introduced;  these  are  defined  in 
28 

such  a  manner  as  to  automatically  satisfy  the  continuity  conditions 
implicit  in  Equations  (86)  and  (87),  viz. 


u  (r,0)  =  -(l/r2sin0)  W/30)  ;  uQ(r,0)  =  (l/rsin0)  (3<J>/3r)  ,  (92) 

u*(r , 0)  =  -(l/r2sin0) (3^*/30)  ;  u*(r,0)  =  (l/rsin0) ( 3^*/3r)  .  (93) 

r  t) 

On  substituting  these  expressions  into  the  expanded  forms  of  Equations 
(88)  and  (89)  the  following  partial  differential  equations  can  eventually 


be  obtained: 


.  ' 
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annular  region:  E2(E2i|j)  =0  R  <  r  <  S  ,  (94) 

exterior  region:  ~(1/k)E2^*  +  E2(E2\Jj*)  =0  S  <  r  <  00  .  (95) 

28  42 

In  these  Equations  E2  denotes  the  Spherical  Harmonic  Operator  *  , 

defined  by: 

E2  e  ( 32/3r2)  +  (1/r2) (32/302)  -  (co t0/r2) ( 3/30)  .  (96) 


In  order  to  describe  the  hydrodynamic  conditions  prevailing 
throughout  the  entire  domain  of  the  modelled  system  it  will  be 
necessary  to  determine  respective  solutions  of  Equations  (94)  and 
(95)  which  satisfy  physically  rational  boundary  conditions. 

II. 4. 2  STIPULATED  BOUNDARY  CONDITIONS 

Recognizing  that  the  reference  sphere  is  impermeable,  and  that 
there  is  no  slip  over  the  surface  thereof,  implies  that: 


ur(R+,8)  =  0  ,  (97) 

u0(R+,6)  =  0  .  (98) 


From  considerations  of  equilibrium  and  continuity  at  the  interface 
separating  the  annular  and  exterior  regions  it  is  necessary  that  the 
pressure,  the  tangential  shear  stress  and  the  velocity  distributions 
within  these  regions  match  identically  in  the  limit,  viz. 


P(s",e)  =  p*(s+,e)  , 

(99) 

t(s",0)  =  x*(s+,e)  , 

(100) 

ur(s",0)  =  u*(S+,0)  , 

(101) 

ue(S",6)  =  u*(S+,0)  . 

(102) 

' 
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Moreover,  the  velocity  vector  ju*  at  any  station  sufficiently  far 
removed  from  the  reference  sphere  must  approach  that  of  the  mainstream, 
U* ,  thus : 

Limit  _u*(r,0)  =  U*  ,  (103) 

r  -*  00 

where  U*  =  [U*,0,0]  in  Cartesian  coordinates  [x,y,z]. 
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II. 5  SOLUTION  OF  THE  DEFINING  EQUATIONS 

II. 5.1  THE  GLOBAL  STREAMFUNCTION  DISTRIBUTION 

The  flowfield  throughout  the  modelled  system  will  be 
completely  defined  by  the  respective  solutions  of  Equations  (94)  and 
(95)  which  satisfy  the  boundary  conditions  stipulated  in  Equations 
(97)-(103).  As  detailed  in  Appendix  C  these  solutions  are  found  to  be: 

iKx»6)  =  (kU*/2)[A/x  +  Bx  +  Cx2  +  Dx4]sin2e  a  <  x  <  3  ,  (104) 

^*(X,0)  =  (kU*/2)[e/X  +  X2  +  Ge"X(l+x)/x]sin29  3  <  X  <  »  (105) 

where  x  denotes  the  normalized  radial  coordinate  defined  by: 

X  =  r//ic  .  (106) 

The  parameters  a  and  3  denote  the  specific  values  of  x  defined  by: 

a  =  R/  /k  ,  (107) 

3  =  S//k  ,  (108) 

and  A,B,C,D SE  and  G  are  expressed  in  terms  of  a  and  3  as  follows: 

A  =  6a3(-2B6-735-15B4-1533+3B4a2-B3a3+333a2+32a3)/J(a,B)  ,  (109) 

B  =  6a (636+2135+4534+453 3-534a2-533ot2-3a5-a5) /J (a, 3)  ,  (110) 

C  =  3(-836-2835-6034-6032+533a3-532a3+33a5+3a5)/J(a,3)  ,  (111) 

D  =  3(484+4B3-383cH-3B2a-Ba3-a3)/J(a,8)  ,  (112) 

E  =  2(-489-2438-60B7-6036+938a+4587a-1036a3+12686a-3035a3 

+21635a+984a5-6034a3+27031+a-433a6+933a5 

-6033a3+27033a-68a6-6ae)/J(a,8)  ,  (113) 

G  =  6el3(2036-2735a+583a3-9033a+2a6)/J(a,3)  ,  (114) 


• 

. 
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where : 

J(oi,,3)  =  (  — 43^  ~ 243  3-180  3  ^  -1803  3+93  3  ot+453  ^a— 10  3  3ct  3 

+18033a-3032a3+93a5-4ct6+9a5)  .  (115) 

The  parameters  a  and  3  appearing  above  are  related  according  to 
Equations  (107),  (108)  and  (83),  thus: 

3  =  a(S/R)  =  a(l-e)-1/3  .  (116) 

II. 5. 2  DISTURBANCE  INTRODUCED  BY  THE  MODELLING  PROCEDURE 

Typical  streamlines,  as  computed  from  Equations  (104)  and 
(105)  ,  are  displayed  in  Figure  17  for  the  modelled  system  possess¬ 
ing  the  representative  porosity  e  =  0.7  (for  which  S/R  =  1.494 
from  Equation  (83)  ,  a  =  3.708  from  Table  3  to  follow,  whence 
3  =  5.540  from  Equation  (116)).  The  streamlines  for  other  porosities 
are  very  similar  in  appearance. 

The  most  important  implication  of  Figure  17  is  that  the 
disturbance  of  the  mainstream  flowfield  introduced  by  the  modelling 
procedure  diminishes  rapidly  with  increasing  distance  from  the  unit 
cell:  this  disturbance  is  effectively  confined  to  a  region  concentric 
with  the  unit  cell  and  possessing  twice  its  radius  (it  will  be 
recalled  that  for  the  corresponding  case  of  diffusion  (Section  1.3.4) 
this  disturbance  was  wholly  confined  to  within  the  unit  cell) .  This 
observation  counters  any  argument  that  the  modelling  procedure  might 
cause  anything  more  than  an  (inevitable)  localized  disturbance  of  the 


mainstream  flowfield. 


.  ' 
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I 


FIGURE  17.  TYPICAL  STREAMLINES  FOR  INCOMPRESSIBLE  CREEPING  FLOW  THROUGH  THE 
MODELLED  SYSTEM  (For  The  Representative  Porosity  e  =  0.7) 
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II. 5. 3  THE  RESISTANCE  OFFEREE  BY  THE  REFERENCE  SPHERE 

Using  the  derived  expression  for  the  s treamf unction  within  the 
annular  region,  i|j(x»0)  ,  it  is  straightforward,  although  tedious,  to 
develop  expressions  for  the  pressure  and  shear  stress  distributions 
over  the  surface  of  the  reference  sphere,  p(a,0)  and  x(a,0)  respec¬ 
tively.  The  normal  and  tangential  resistance  forces  offered  by  the 
reference  sphere  (F  and  F^  respectively)  may  then  be  calculated  by 
integrating  these  distributions  over  the  entire  surface  thereof.  The 
total  resistance  force,  F,  offered  by  the  reference  sphere  is  so 
derived  in  Appendix  Cs  viz. 

F  =  F  +  Fq  =  6ttijU*R  £(a,0)  ,  (117) 

r  0 

where  the  function  £(a,0)  is  defined  by: 

£(a,3)  =  4(-636-213 5-450t+-453 3+50l+a2-4  53 3a2+3a5+a5) /J(a,3)  .  (118) 


It  will  be  demonstrated  later  (Table  3)  that  as  e  -*  1.0,  £(a,3)  ->  1.0 

28 

so  that  Equation  (117)  then  approaches  Stokes  Law  for  creeping  flow 
past  a  single  sphere  in  space,  as  would  be  expected,  viz. 


F 


Stokes 


6ttuU*R  . 


(119) 


The  term  £(a,3)  therefore  represents  the  extent  to  which  the  spheres 
surrounding  the  reference  sphere  increase  the  latter's  resistance  over 
that  predicted  by  Stokes  Law  for  a  single  isolated  sphere. 


' 
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II.  5.  4  THE  TOTAL  RESISTANCE  OFFERED  BY  A  SWARM  OF  MONOSIZED  SPHERES 

This  analysis  has  heretofore  been  concerned  with  an  unbounded 
homogeneous  swarm  of  spheres.  As  noted  earlier  in  Section  II. 2  the 
Darcy  Equation  (81)  is  valid  for  macroscopically  uniform  flow  through 
such  a  system.  According  to  this  equation  the  total  resistance  offered 
by  a  swarm  containing  N  monosized  spheres  of  radius  R  will  be  (Equation 
(C71))  : 

F  =  6ttuU*KN  { 2ot2/9  (1-e)  }  .  (120) 

The  total  resistance  offered  by  the  modelled  system  may  be 
evaluated  by  summing  Equation  (117)  over  all  IT  spheres  (since  each 
individual  sphere  may  be  regarded  as  the  reference  sphere  in  turn) , 
thus  : 

F  =  6ttiiU*RN  5(a,0)  .  (121) 

Since  the  modelled  system  is  to  be  quantitatively  representative  of 
the  original  swarm  of  monosized  spheres  then  the  resistance  forces 
predicted  by  Equations  (120)  and  (121)  must  be  identical.  This 
generates  the  condition  that: 

2a2  -  9(l-e)S(a,3)  =  0  ,  (122) 

with  the  parameters  a  and  3  being  related  according  to  Equation  (116) , 
thus  : 

3  =  a(l-e) -1 /3  .  (123) 

Consequently,  Equation  (122)  may  be  re-expressed  as  the  following 
implicit  function  of  a  and  e,  viz. 


'* 
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2a2  -  9 (1-e) £ (a ,a(l-e)  ^/3)  =  0  .  (124) 

This  equation  effectively  defines  the  functional  relationship  existing 
between  a  =  R/ /k  and  the  porosity  e  for  any  homogeneous  swarm  of  mono¬ 
sized  spheres. 

II. 5. 5  COMPUTED  RESULTS  FOR  MONOSIZED  SPHERES 

Recalling  from  Equation  (118)  the  extremely  complicated  nature 
of  the  function  £(a,8)  it  is  apparent  that  a  cannot  be  extracted  from 
Equation  (124)  above  as  a  function  of  e  in  closed  form.  Resort  must 
therefore  be  made  to  an  iterative,  and  in  consequence  pointwise, 
technique  of  solution.  As  expected,  Equation  (124)  constitutes  a 
single  valued  function  over  the  entire  porosity  range;  the  regula  falsi, 
technique  provides  a  rapid  iteration  on  a  using  a  digital  computer. 

R.ep  resent  at  ive  values  of  the  a  =  a(e)  relationship,  satisfying 
Equation  (124)  to  within  better  than  10  ^ ,  were  obtained  using  this 
technique.  These  results  are  presented  (to  four  significant  figures) 
in  Table  3  for  reference,;  it  will  be  noted  that  these  predictions  are 
physically  consistent  at  both  porosity  limits,  thus  as  e  -*  0 ,  a  *>  00 
ana  as  e  -*-1.0,  a->0. 

It  is  particularly  instructive  to  compare  the  resistance  of 
the  original  system  as  predicted  by  Stokes  Law,  when  assuming  each 
sphere  to  be  hydrodynamically  independent  of  the  remainder,  with  that 
predicted  by  the  presented  theory.  This  ratio,  VT,  for  monosized 
spheres  is  given  by  (Equation  (C86)): 


W  =  9(l-e)/2a2  =  l/£(a,a(l-e)  x/3)  . 


(125) 


' 

. 
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TABLE  3 

THE  PREDICTED  DEPENDENCE  OF  a  AND  W 
ON  e  FOR  HOMOGENEOUS  SWARMS  OF 
MONO SIZED  SPHERES 


E 

a 

W  =  1/S(a,3) 

1.0 

0.0 

1.0 

0.9999999999 

0. 2003xl0~4 

0.9994 

0.99999999 

0.2124xl0”3 

0.9968 

0.999999 

0.2137xl0~2 

0.9850 

0.9999 

0 . 2200xl0_1 

0.9300 

0.999 

0 . 7282xl0_1 

0.8487 

0.99 

0.2584 

0.6738 

0.95 

0. 7071 

0.4501 

0.90 

1.185 

0.3207 

0.85 

1.684 

0.2379 

0.80 

2.247 

0.1782 

0.75 

2.908 

0.1331 

0.70 

3.  708 

0.9819xl0_1 

-1 

0.65 

4. 706 

0.7112x10 

-1 

0.60 

5.986 

0.5023x10 

-1 

0.55 

7.677 

0.3436x10 

-1 

0.50 

9.974 

0.2262x10 

-1 

0.45 

13.19 

0.1423x10 

-2 

0.40 

17.83 

0.8496x10 

-2 

0.35 

24.74 

0.4780x10 

-2 

0.30 

35.44 

0.2508x10 

-2 

0.25 

52.93 

0.1205x10 

-3 

0.20 

83.87 

0.5118x10 

-3 

0.15 

146.0 

0.1794x10 

0.10 

302.7 

0.4419xl0"4 

-5 

0.05 

967.9 

0.4564x10 

0.0 

00 

0.0 
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The  W  “  W(e)  relationship  predicted  by  this  equation  is  displayed  in 
Table  3. 

II. 5. 6  COMPARISON  OF  COMPUTED  RESULTS  WITH  EXPERIMENTAL  DATA 

Figure  19  displays  the  a  =  a(e)  relationship  for  monosized 
spheres,  as  predicted  by  Equation  (124),  together  V7ith  representative 

or  /Q  / 7  /Q  cn 

data  from  the  literature  *  *  *  5  .  The  overall  agreement  can 

be  seen  to  be  very  satisfactory. 

In  keeping  with  literature  trends  the  preceding  results  are 

presented  in  alternative  form  in  Figure  20,  which  displays  the 

W  =  W(e)  relationship  predicted  by  Equation  (125) . 

In  connection  with  both  Figures  19  and  20  it  should  be  emphasized 

43  48 

that  the  data  of  Happel  and  Epstein  and  Martin  et  al  was  obtained 

40 

using  regular  assemblages  of  spheres  rather  than  using  a  random 

37 

assemblage  thereof  .  This  is  possibly  the  reason  as  to  why  their 
data  points  lie  rather  less  close  to  the  present  predictions  than  do 
those  of  the  other  cited  workers  (who  confined  their  attention  to 
randomly  arranged  spheres) . 

Present  Eooperimental  Work 

In  order  to  confirm  the  validity  of  existing  experimental  data 
in  the  lower  porosity  range  a  number  of  permeability  determinations 
were  carried  out  using  a  cylindrical  brass  cell  which  was  filled, 
under  distilled  water,  with  monosized  stainless  steel  spheres  (Figure 
18)  .  Distilled  water  at  a  constant  temperature  could  be  made  to  flow 
through  the  cell  at  a  pre-selected  constant  flowrate  (effected  by 
means  of  a  commercial  constant-flow  regulator).  For  each  constant 


(  ) 


. 
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flowrate  the  pressure  gradient  in  the  direction  of  flow  was  ascertained 

by  means  of  a  pressure  transducer  (periodically  calibrated  against  a 

high  precision  differential  manometer) .  The  permeability  of  the  system 

was  then  calculated  according  to  Equation  (81) .  Finally,  the  porosity 

of  the  system  was  determined  by  a  weighing  technique  (although  the 

minimum  porosity  obtainable  with  monosized  spheres  is  0.2575, 

40 

corresponding  to  a  regular  rhombohedral  array  ,  it  is  seldom  pos¬ 
sible  to  obtain  a  porosity  lower  than  0.37  with  randomly  arranged 
spheres) . 

Several  measurements  were  carried  out  within  the  Reynolds 
number  range  Re  =  0.08  -  1.40  (based  on  the  superficial  liquid 
velocity  and  the  particle  diameter)  .  Within  this  range  little 
discernible  change  in  the  measured  permeability  could  be  observed. 

This  data  is  presented  in  normalized  form  in  Table  4  below  and  in 
Figure  23  to  follow.  As  displayed  in  Figures  19  and  20  this  data 
conforms  closely  with  previously  reported  experimental  data. 


f 


■  - 


TABLE  4 


EXPERIMENTAL  PERMEABILITY  DATA  FOR  AN 
HOMOGENEOUS  SWARM  OF  MONOSIZED  SPHERES 


e  =  0.379 

Re 

a 

W 

0.081 

19  .08 

0.00768 

0.093 

19.12 

0.00765 

0.168 

19 .10 

0.00766 

0.192 

19.09 

0.00769 

0.251 

19.12 

0.00765 

0.290 

19.10 

0.00766 

0.552 

19.10 

0.00766 

0.670 

19.13 

0.00764 

0.807 

19.11 

0.00765 

0.932 

19 .14 

0.00763 

1.38 

19 .16 

0.00762 
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C  CYLINDRICAL  BRASS  CELL 

S  STAINLESS  STEEL  SPHERES  (1mm  die.  bearing  balls) 

Q  DISTILLED  WATER  SUPPLY  (  20 -100  psig.) 

P  PRESSURE  TAPS  (spaced  3"  apart) 

W  WEIGHING  SCALE 

R  CONSTANT  FLOW  REGULATOR 
(Moore  Instrument  Co.  model  63  BU) 

T  PRESSURE  TRANSDUCER  (Whittaker  Corp.  model  KP15) 

D  DIGITAL  VOLTMETER  (Non-Linear  Systems  Inc.  model  X-3) 

FIGURE  18.  EQUIPMENT  USED  FOR  THE  EXPERIMENTAL  DETERMINATION 
OF  THE  PERMEABILITY  OF  A  PACK  OF  SPHERES 
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a 


FIGURE  19. 


the  ?rit: cm  dipestesxz  ot  ^  oh  £  for  >*:nosized 

5 p HIRES :  COMPARISON  -171  EXPERIMENTAL  DATA 
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II. 6  DISCUSSION  OF  PREVIOUS  WORK 

Numerous  experimental  and  theoretical  studies  concerning  the 

present  problem  have  been  reported  in  the  literature.  As  these 

52 

have  been  reviewed  in  conscientious  detail  by  Scheidegger  and  Happel 
42 

and  Brenner  discussion  will  here  be  limited  to  those  theoretical 
treatments  which  are  considered  to  be  of  principal  interest  and 
significance. 

The  first  worker  to  establish  a  model  concerning  the  problem 

36 

of  interest  appears  to  have  been  Cunningham  who,  in  1910,  studied  the 
rate  of  sedimentation  of  a  cloud  of  spheres  in  a  fluid  medium.  He 
postulated  that  each  particle  within  the  cloud  would  effectively  be 
limited  to  motion  within  a  concentric  mass  of  fluid.  He  further 
assumed  the  outer  boundary  of  this  fluid  envelope  to  be  solid,  presum¬ 
ably  corresponding  to  the  surface  of  the  surrounding  spheres.  Although 
of  fundamental  significance,  this  model  presents  difficulty  in  that 
the  dimensions  of  the  spherical  envelope  must  be  fixed  by  additional 
empirical  considerations. 

Based  on  Kozeny’s^5  classic  capillary  model.  Carman3^  developed 
the  following  semi— empirical  correlation  for  the  permeability,  k,  of 
a  bed  of  monosized  spheres : 

k  =  R2e3/(9K  (1-e)2}  ,  (126) 

c 

which  together  with  Equations  (107)  and  (125)  becomes. 

W  =  e 3/{ 2K  (l-e)>  • 
c 


(127) 


, 
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The  parameter  denotes  the  so  called  'Carman  constant’  which  has  to 
be  determined  by  experiment.  The  universally  acknowledged  value  of 
for  spheres  is  5.0,  for  which  Equation  (127)  becomes  reasonably 
representative  of  a  vast  quantity  of  experimental  data  within  the 
porosity  range  0.26  <  e  <  0.7  (Figure  20). 

One  of  the  more  widely  quoted  theories  concerning  the  present 

29 

problem  is  that  due  to  Brinkman  .  He  extended  the  model  proposed  by 
33 

Bruggemann  for  electrical  conduction  (Section  1.4)  to  the  case  of 
creeping  flow,  the  representation  underlying  both  treatments  being 
that  of  a  spherical  particle  embedded  within  a  uniform  porous  mass. 
Brinkman  obtained  the  following  solution  for  an  homogeneous  swarm  of 
monosized  spheres; 

W  =  1  +  a  +  a2/3  =  1  -  (3/4) (1-e) [/8/(l-e)-3  -l]  .  (128) 

However,  this  solution  is  entirely  unsatisfactory  for  e  <0.4  and,  in 
fact,  predicts  zero  permeability  for  e  =  1/3  (Figure  20) t. 


t  Preferring  here  to  the  predictions  of  the  present  theory 
it  is  elucidating  to  consider  the  limit  of  Equation  (125)  as 
3  ->  a,  that  is  as  the  thickness  of  the  annular  region  in  the  proposed 
model  tends  to  zero.  Under  these  conditions  the  proposed  model 
reduces  to  the  simpler  one  employed  by  Brinkman;  an  inspection  of 
Equation  (125)  in  conjunction  with  Equation  (118)  ultimately  yields; 

Limit  W  =  Limit  [l/6(a,3)]  =  1  +  a  +  a2/3  . 

3  ->  a  3  -*  a 

This  is  identical  with  the  Brinkman  solution  as  would  be  expected. 

Aiming  to  generalize  his  model  to  include  porous  spheres 
Brinkman^  made  attempts  to  incorporate  an  annular  region  which 
could  be  varied  in  size  so  as  to  present  the  experimental  facts 
throughout  the  entire  porosity  range.  However,  this  valuable  contri¬ 
bution  has  been  neglected  in  the  literature  to  date,  probably  on 
account  of  the  fact  that  his  revised  solution  was  physically 
inconsistent  and  was  presented  without  derivation. 


' 

* 
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, 


II. 6] 


84 


56 

Uchida  attempted  a  rigorous  solution  of  the  Kavier-S tokes 

and  continuity  equations  for  incompressible  creeping  flow  through 

an  unbounded  cubic  assemblage  of  monosized  spheres.  However,  the 

problem  of  simultaneously  satisfying  the  appropriate  boundary  conditions 

on  spherical  as  well  as  cubical  surfaces  proved  to  be  intractable  and 

he  was  only  able  to  extract  an  approximate  solution;  this  failure  to 

solve  exactly  the  specified  boundary  value  problem  resulted  in  very 

unsatisfactory  predictions. 

41 

Happel  proposed  an  interesting  cell-type  model  which  is  now 

referred  to  as  the  ’free  surface  model’.  It  is  apparently  based  on  a 

3  6 

refinement  of  Cunningham’s  previously  discussed  cell  model,  the  solid 
outer  envelope  thereof  here  being  replaced  by  a  free  fluid  surface  on 
which  both  the  normal  velocity  and  the  tangential  shear  stress  vanish. 
Using  this  representation  Happel  derived  the  following  expression  for 
monosized  spheres: 

W  =  [6  -  9(1- e)1/3  +  9(l-e)5/3  -  6(l-e)2]/[6  +  4(l-e)5/3]  .  (129) 

Although  the  free  surface  model  may  be  superficially  unconvincing  the 

agreement  of  its  predictions  with  experimental  data  is  undeniable. 

It  is  interesting  to  note  that  the  predictions  of  the  presented 

theory  are  in  remarkably  close  agreement  with  those  of  Happel  for 

e  >  0.7,  nowhere  differing  by  more  than  1%  (Figure  20). 

A6 

Kuwabara  modified  the  free  surface  model  by  replacing  the 
zero  shear  stress  condition  at  the  outer  envelope  by  one  of  zero 
vorticity.  However,  his  predictions  are  far  less  satisfactory  than 
those  of  Happel. 
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Various  statistical  models  have  been  advanced  to  study  fluid 

flow  through  porous  media  in  general.  A  number  of  these  have  been 

52 

reviewed  by  Scheidegger  ^  (who  himself  proposed  a  model  in  which  the 
theory  of  Brownian  motion  is  applied  to  the  problem  in  question) . 

Of  significant  interest  amongst  these  is  the  approach  of  Yuhara  who 
postulated  an  analogy  between  laminar  flow  in  porous  media  and 
turbulent  flow  in  bulk  fluids.  On  a  somewhat  different  theme  Broadbent 
and  Hammersley  proposed  a  model  in  which  any  ’randomness’  is  attached 
to  the  medium  rather  than  to  the  fluid  (this  is  generally  the  case  in 
practice) .  Although  such  approaches  are  extremely  interesting  they 
tend  (at  present)  to  be  too  qualitative  in  nature  to  permit  any  direct 
engineering  application. 


’ 


II. 6] 


86 


6 


FIGURE  20. 


THE  DEPENDENCE  OF  W  ON  e  FOR  MONOSIZED  SPHERES: 
COMPARISON  WITH  PREVIOUS  WORK 
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IX. 7  THE  EFFECTS  OF  THE  SIZE  DISTRIBUTION 

It  was  noted  in  Section  1.3  that  for  diffusive  flow  processes 
occurring  within  an  homogeneous  swarm  of  spheres  the  size  distribution 
had  little  or  no  effect  on  the  diffusivity  factor.  In  contrast,  for 
fluid  flow  the  presented  theory  predicts  a  definite  dependence  of  the 
permeability  on  the  prevailing  size  distribution  of  the  spheres  (Appen¬ 
dix  C) .  The  specific  predictions  which  follow  have  been  confined  to 
binary  mixtures  of  spheres  in  order  to  permit  a  more  meaningful  inter¬ 
pretation  of  the  effects  of  the  size  distribution. 

Size  Ratio  5:1 

Figure  21  displays  the  W  =  W(e)  relationships  predicted  by 
Equations  (C84)  and  (C82)  for  binary  mixtures  of  spheres  possessing 
arbitrary  radii  R  and  R/5,  in  the  relative  proportions  100%,  50%,  25%, 
10%,  1.74 %  and  0%  by  numbers  of  the  larger  species.  As  would  be  expec¬ 
ted  the  results  for  100%  and  0%  are  identical  since  both  systems  are 
then  composed  of  monosized  spheres.  Mixtures  containing  more  than  50% 
of  the  larger  species  exhibit  characteristics  very  close  to  those  for 
monosized  spheres  and  have  therefore  been  omitted  from  Figure  21  to 
preserve  clarity. 

The  solutions  for  50%,  25%  and  10%  display  an  increasing  depar¬ 
ture  from  that  for  monosized  spheres  and,  in  fact,  a  maximum  departure 
is  predicted  for  1.74%.  As  this  proportion  further  decreases  to  0% 
the  solution  rapidly  re-approaches  that  for  monosized  spheres.  Hence, 
the  solution  for  any  arbitrary  mixture  possessing  a  5:1  size  ratio  will 
always  lie  between  that  for  monosized  spheres  and  that  for  the  mixture 
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containing  1.74%  of  the  larger  species. 

Size  Ratio  2:1 

Figure  22  displays  the  results  predicted  by  Equations  ( C 8 4 ) 
and  (C82)  for  binary  mixtures  of  spheres  possessing  arbitrary  radii 
R  and  R/2.  However,  only  the  two  limiting  solutions  (viz.  that  for 
monosized  spheres  and  that  exhibiting  the  maximum  departure  therefrom) 
have  been  included  because,  although  the  general  trends  are  closely 
similar  to  those  for  the  5:1  size  ratio,  the  corresponding  departures 
are  very  much  less  pronounced.  The  maximum  departure  is  here  exhibited 
by  the  mixture  containing  15.2%  of  the  larger  species. 

The  most  important  inference  to  be  drawn  from  Figures  21  and 
22  is  that  W,  and  hence  the  permeability  of  the  system,  do  indeed 
depend  on  the  prevailing  size  distribution  of  the  spheres;  this 
dependency  increases  rapidly  with  the  size  ratio.  Equations  (C84) 
and  (C82)  may  likewise  be  employed  to  accomodate  ternary,  quaternary 
and  higher  order  mixtures  should  the  occasion  arise. 
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FIGURE  21.  THE  EFFECTS  OF  THE  SIZE  DISTRIBUTION  IN  BINARY 
MIXTURES  OF  SPHERES  (SIZE  RATIO  5:1) 


II. 7] 


90 


€ 


FIGURE  22.  THE  EFFECTS  OF  THE  SIZE  DISTRIBUTION  IN 

BINARY  MIXTURES  OF  SPHERES  (SIZE  RATIO  2:1) 
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II. 8  THE  EFFECTS  OF  THE  REYNOLDS  NUMBER 

The  preceding  theoretical  results  have  all  been  based  on  the 
premise  of  creeping  flow,  viz.  that  inertial  forces  can  be  neglected 
in  comparison  with  viscous  forces.  In  practice  this  is  generally  a 
valid  assumption  for  liquid  flow;  however,  exceptions  may  be  encountered 
from  time  to  time.  To  this  end,  many  investigations  have  been  directed 
towards  determining  a  universal  Reynolds  number,  above  which  the  creep¬ 
ing  flow  assumption  ceases  to  be  strictly  valid  for  liquid  flow  through 
porous  media  in  general.  For  swarms  of  monosized  spheres,  in  particular, 
this  specific  Reynolds  number  is  generally  acknowledged  to  be  about 

unity  (Figure  23) ,  when  based  on  the  superficial  velocity  of  the  fluid 

32  38  49  51 

and  the  particle  diameter  *  *  *  However,  no  universal  result 

concerning  porous  media  in  general  has  so  far  been  forthcoming. 

It  is  possible  that  this  specific  number  for  an  homogeneous 
swarm  of  spheres  could  be  determined  by  means  of  the  proposed  model. 

This  has  been  attempted  by  accounting  for  the  inertial  force  term, 
pu.Vu,  in  the  Navier-Stokes  Equation  (84),  which  describes  conditions 
within  the  annular  region,  and  by  introducing  an  equivalent  macroscopic 
term,  pu*.Vu*,  into  the  Brinkman  Equation  (85),  which  describes 
conditions  within  the  exterior  region,  thus: 

annular  region:  yV2u  =  Vp  +  pu.Vu  R  <  r  <  S  ,  (130) 

exterior  region:  -(p/k)u *  +  pV2u*  =  Vp*  +  pu*.Vu*  S  <  r  <  00  .  (131) 

The  above  equations  have  defied  analytical  attempts  of  solution.  Perturb- 
ation  techniques  or  numerical  methods  would  appear  to  be  inevitable. 
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FIGURE  23.  THE  DEPENDENCE  OF  a  ON  THE  REYNOLDS  NUMBER  FOR  MONOSIZED  SPHERES:  EXPERIMENTAL  DATA 
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II. 9  THE  EFFECTS  OF  PARTICLE  POROSITY 


As  noted  previously  (Section  1.3.6)  porous  particles  frequently 
play  an  important  role  in  catalyst  beds.  The  problem  often  arises  as 
to  the  fraction  of  a  flowing  fluid  which  passes  through  a  single 
particle  contained  witr.ir.  such  a  bed.  With  this  in  mind  it  is 
particularly  instructive  to  calculate  the  streanf unction  distribution 
for  incompressible  creeping  flow  through  and  around  a  single  isolated 
porous  sphere  (of  radius  R  and  permeability  <) .  Incidentally,  this 
analysis  will  serve  to  demonstrate  the  effectiveness  of  the  Brinkman 
Equation  (82)  in  solving  an  hitherto  intractable,  two-region  problem. 

II.  9.1  FUNDAMENTAL  DIFFEREHTIAL  EQUAIIONS 

By  analogy  with  Equations  (94)  and  (95) ,  the  equations  describing 
the  flovfield  within  the  porous  sphere  and  within  the  surrounding  fluid 
region  will  be : 

porous  sphere:  -(l/<)E,dt*  +  E‘-(E"t*)  =0  0  <  r  <  R  ,  (132) 

surrounding  fluid:  E-(E^v)  =0  R  <  r  <  «=  .  (133) 


II.  9. 2  STIPULATED  BOUND API  I0EIITI0DS 


The  following  boundary  conditions  must  here  be  stipulated,  viz. 


p*(R  ,6)  =  p(R  ,9)  , 

(134) 

-r*(R~,e)  =  t(R+,9)  , 

(135) 

uJ(R  ,9)  =  ur(R+ ,e)  , 

(136) 

u$(R  ,9)  “  u4(R+ ,9)  , 

(137) 

Limit  ti(r ,  0)  =  _U  , 

(138) 

' 
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where  tJ  ==  [U,0,0]  denotes  the  velocity  vector  of  the  undisturbed  main¬ 
stream;  the  arguments  on  which  the  above  boundary  conditions  are  based 
are  presented  in  Section  II. 4. 2.  In  addition,  a  restriction  demanding 
finiteness  of  the  flowfield  must  be  imposed,  in  particular  that 


u£(O,0)  must  be  finite.  (139) 

II. 9. 3  SOLUTION  OF  THE  DEFINING  EQUATIONS 


The  solutions  of  Equations  (132)  and  (133)  which  satisfy  the 
boundary  conditions  implicit  in  Equations  (134) -(139)  can  be  shown  to 
be : 

^*(X>6)  =  (kU/2)[Mx2  +  N(sinhx-Xcoshx)  /x]sin20  0  <  x  <  v  >  (140) 

lKx»0)  =  (kU/2)  [A/x  +  Bx  +  x2 ] sin20  v  <  X  <  00  »  (141) 


whe  re : 

X  =  r/ ,  (142) 

v  =  R//k"  ,  (143) 

and : 

M  =  (tanhv-v) /(tanhv-v-2v3/3)  ,  (144) 

N  =  (2v3/coshv) / (tanhv-v-2v3/3)  ,  (145) 

A  =  2v3 (tanhv{l+v2/2}-v{l+v2/6}) /  (tanhv-v-2v3 /3)  ,  (146) 

B  =  -v3 (tanhv-v) /( tanhv-v-2v3/3)  .  (147) 


When  the  sphere  ceases  to  exist  (R  =  0,  v  =  0),  Equation  (141)  reduces 

to  the  limiting  form  \p  =  (l/2)Ur2sin20 .  This  is  the  familiar  stream- 

2  8 

function  representation  of  an  undisturbed  flowfield,  U,  as  would  be 


expected. 
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II. 9. 4  THE  RESISTANCE  OFFERED  BY  AH  ISOLATED  POROUS  SPHERE 

The  resistance  force,  F,  offered  by  an  isolated  porous  sphere 
may  be  calculated  by  means  of  Equation  (C62)  ,  which  here  assumes  the 
simpler  form: 


F  =  6ttpURC12A  -  2Bv2)  /9v 3  .  (148) 

Substituting  for  A  and  B  from  Equations  (146)  and  (147)  ultimately 
yields  A  Generalized  Form  of  Stokes  Law3  for  a  Porous  Sphere ^  thus : 

F  =  6uyUR  C(v)  ,  (149) 

where : 

1  -  (7/3v)tanhv  +  (4/v2)  -  (4/v3)tanhv 

C(v)  =  -  (150) 

1  +  (3/2v2)  -  (3/2v3)tanhv 

For  solid  spheres  v  =  whence  C(v)  =  1.0;  Equation  (149)  then  becomes 
identical  with  Stokes  Law  as  would  be  expected. 

For  v  >  10,  C(v)  ~  l-(7/3v).  In  general  v  >>  100  for  commercial 

porous  spheres,  whence  C(v)  -  1.0  (Figure  24).  The  immediate  inference 
is  that  for  v  >>  100  an  insignificant  quantity  of  fluid  actually  passes 
through  the  particle  itself;  in  other  words  the  fluid  prefers  to  flow 
around  the  porous  sphere  just  as  if  it  were  solid  (Figure  25).  The 
implication  of  this  for  flow  within  a  swarm  of  porous  spheres  is  that 
most  of  the  fluid  will  flow  between  the  spheres;  only  a  small  fraction 
can  be  anticipated  to  flow  through  the  individual  spheres.  In  view  of 
this  a  bed  of  catalyst  pellets  can  be  expected  to  function  more 
efficiently  when  composed  of  a  multitude  of  very  small  solid  particles 
than  when  composed  of  a  few  much  larger  porous  particles  possessing 
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the  same  overall  volume  (since  catalysis  is  predominantly  a  surface 
phenomenon) . 

The  proposed  model  for  an  homogeneous  swarm  of  porous  spheres 
(Figure  4)  may,  in  theory,  be  extended  to  the  present  case  of  incom¬ 
pressible  creeping  flow.  Unfortunately  the  algebra  involved  grows 
so  excessively  complicated  that  all  attempts  of  completing  the 
solution  have  failed.  However,  in  view  of  the  inferences  drawn  from 
the  above  analysis  for  an  isolated  porous  sphere,  such  a  solution  would 
appear  to  have  limited  applicability. 

II. 9. 5  TEE  RESISTANCE  OFFERED  BY  A  DENSE  CLUSTER  OF  SPHERES 

It  should  be  stressed  that  the  results  developed  above  also  apply 
to  the  case  of  a  single  porous  sphere  sedimenting  within  an  unbounded 
extent  of  fluid.  Of  more  interest,  however,  are  the  sedimentation 
characteristics  of  a  dense  cluster  of  spheres,  particularly  when  the 
cluster  itself  approximates  a  porous  spheret. 

At  present  there  does  not  appear  to  be  available  any  theoretical 

result  relating  to  the  sedimentation  of  such  a  cluster  of  spheres. 

However,  the  results  developed  above  for  a  single  porous  sphere 

may  be  applied  directly  to  the  case  of  a  sedimenting  spherical  cluster 

(of  radius  R  and  permeability  <)  composed  of  monosized  spheres  (of 
c 

t  Whilst  studying  the  behaviour  of  biological  macro-molecules 
in  solution,  Brinkman30  noted  that  such  particles  showed  a  tendency 
to  sediment  in  dense  clusters  which  approximated  porous  spheres.  (He 
attempted31  to  incorporate  this  observation  into  his  previously 
discussed  model  for  an  homogeneous  swarm  of  solid  spheres29  but  was 
unable  to  obtain  a  physically  consistent  solution) . 
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radius  R  ).  The  value  of  v  for  such  a  cluster  will  be: 

o 


V  -  R  /Sk  =  (R  /Jk)  (R  /R  ) 
c  p  c  p 

=  a(R  /R  )  , 
c  p 


(151) 


with  a  being  directly  obtainable  from  Figure  (19),  or  Table  3,  if 

the  porosity  of  the  cluster  is  known  (this  will  generally  be  between 

0.40  and  0.44  for  monosized  spheres  in  contact).  Hence,  a  knowledge  of 

(R  /R  )  and  a  for  the  cluster  is  sufficient  to  evaluate  its  sedimenta- 
c  p 

tion  characteristics  defined  by  Equations  (149)  and  (150). 

Inherent  in  the  above  discussion  of  an  isolated  cluster  is 


the  presumption  that  the  cluster  contains  a  sufficient  number  of 

spheres  to  approximate  a  sphere  itself.  For  monosized  spheres  in 

contact  this  necessitates  that  (R  /R  )  >5  (from  geometric  consider- 

c  p 

ations) . 


*• 


II. 9] 


99 


re 

o 

Z 

o 

od 

re 

H 

§ 

£ 

o 

5ZS 

M 

Pu 

W 

a 

o 

w 

hJ 

PQ 


co 

to 


o 

o 

z 


o 

CM 

II 

r> 

<u 

CO 

<TJ 

o 

(U 

> 

•H 

4J 

CCJ 

•u 

c 

<u 

CO 

(U 

CL 

<u 

Pi 

<D 

JC 

H 

>1 

o 

Pm 


Z 

O 

Pm 

CO 

W 

Z 


H 

CO 


3 

a 

M 

p-i 

c 


a 

w 

re 

pL, 

co 

co 


Z 

O 

P4 

o 


PM 

Q 

w 

H 

< 

►J 

O 


co 


LO 

CM 


w 

§ 


U-t 


I 


100 


II. 10  THE  EFFECTS  OF  PARTICLE  SHAPE  AND  ORIENTATION 

The  effects  of  particle  shape  and  orientation  were  considered 
in  Part  IB  during  the  study  of  diffusion  through  certain  arrangements 
of  oblate  and  prolate  spheroids.  A  similar  analysis  for  incompressible 
creeping  flow  has  been  attempted.  In  essence,  this  problem  consti¬ 
tutes  one  of  obtaining  general  solutions  of  Equations  (94)  and  (95) 
in  spheroidal  coordinates.  Unfortunately,  all  attempts  at  determining 

an  analytic  solution  of  this  latter  equation  have  failed  (although 

A  2 

such  a  solution  of  the  former  equation  -is  available  )  .  A  numerical 
approach,  as  previously  suggested  in  Section  II. 8  for  the  evaluation 
of  inertial  effects,  therefore  appears  inevitable.  However,  it 
seems  likely  that  the  solution  so  obtained  would  suggest  a  perme¬ 
ability  dependence  on  the  eccentricity  and  mode  of  orientation  of 
the  spheroids  similar  to  the  diffusivity  dependence  on  these  parameters 
discussed  in  Part  IB. 
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II .  11  SUMMARY 

The  presented  results  demonstrate  that  the  proposed  model 
offers  a  satisfactory  representation  of,  and  provides  valuable 
physical  insight  into,  incompressible  creeping  flow  occurring 
within  an  homogeneous  swarm  of  spheres.  For  solid  monosized 
spheres  of  known  radius,  the  permeability  may  be  evaluated  by 
iteration  of  Equation  (124)  or  by  interpolation  of  Table  3.  These 
predictions  agree  well  with  experimental  data  over  the  entire 
porosity  range,  lending  heavy  support  to  the  realistic  nature  of 
the  proposed  model  and  to  the  acceptability  of  the  assumptions 
implicit  therein. 

This  analysis  has  demonstrated  that,  for  spheres,  the 
permeability  is  not  invariant  with  the  size  distribution,  inferring 
that  fluid  flow  measurements  can  be  expected  to  yield  quantitative 
information  relating  to  pore  sizes,  pore  size  distributions  and 
specific  surface  areas  of  porous  media  in  general;  these  conclusions 
are  also  in  accord  with  experimental  observations. 

Although  the  presented  results  have  been  derived  specifically 
for  creeping  flow,  there  exists  a  great  deal  of  experimental  data 
purporting  to  their  validity  up  to  a  Reynolds  number  of  at  least 
unity.  Since  this  is  substantially  higher  than  that  usually 
encountered  with  liquid  flow  in  porous  media,  the  preceding  results 
may,  in  such  cases,  be  applied  with  confidence. 
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CONCLUSION 

The  nature  of  the  results  derived  in  this  work  serve  to  demon¬ 
strate  that  the  proposed  model  satisfactorily  fulfils  the  original 
objectives  of  (1)  being  able  to  predict  successfully  diffusive  as 
well  as  hydrodynamic  flow  processes  when  occurring  individually  with¬ 
in  an  homogeneous  swarm  of  spheres,  and  (2)  providing  valuable 
physical  insight  into  these  processes  when  occurring  within  uncon¬ 
solidated  porous  media  in  general. 

The  diffusivity  (conductivity)  factor  of  an  homogeneous  swarm 
of  spheres  may  be  evaluated  if  the  porosity  is  known.  When  the  size 
distribution  is  known  the  permeability  also  may  be  evaluated.  These 
predictions  display  similar  trends  to  experimental  data  over  the 
entire  porosity  range,  observations  which  justify  having  resorted  to 
the  same  mathematical  model  throughout  this  work. 

The  diffusivity  (conductivity)  factor  has  been  demonstrated  to 
be  invariant  with  the  size  distribution  of  the  spheres,  thereby  implying 
that  such  measurements  cannot  be  expected  to  provide  quantitative 
information  relating  to  pore  sizes,  pore  size  distributions  or 
specific  surface  areas  of  porous  media  in  general.  In  contrast, 
however,  the  permeability  has  been  demonstrated  to  be  significantly 
dependent  on  the  prevailing  size  distribution.  Hence,  caution  is 
advisable  when  comparing  theoretical  or  experimental  permeability  data, 
obtained  for  a  given  system,  with  experimental  data  obtained  from 
systems  having  like  porosities  but  different  or  unknown  size 


distributions . 
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The  important  effects  of  particle  shape  and  orientation  were 
assessed  during  the  study  of  diffusion  (conduction)  through  systems 
of  co-axially  orientated  oblate  and  prolate  spheroids.  The  proposed 
model  also  yielded  encouraging  predictions  for  these  fundamentally 
important  anisotropic  systems. 

To  summarize  then,  the  proposed  model  provides  valuable 
insight  into  the  understanding  of  diffusion  (conduction)  as  well 
as  fluid  flow  when  occurring  individually  within  unconsolidated 
porous  media.  It  is  totally  different  in  concept  from  the  classical 
capillary  model  which,  although  widely  utilized,  remains  rather 
unrealistic  and  restrictive  in  practice  (its  most  serious  limitation 
being  that  it  is  wholly  unproductive  in  the  study  of  diffusive  flow 
processes);  indeed,  most  unconsolidated  porous  media  (and  many 
consolidated  ones)  are  more  closely  approximated  by  an  homogeneous 
swarm  of  spheres  or  co-axially  orientated  spheroids  than  by  a  bundle 
of  non- interconnected  capillaries . 

With  reference  to  further  work  it  is  suggested  that  the  proposed 
model  be  extended  to  the  study  of  diffusive  and  hydrodynamic  flow 
processes  when  occurring  simultaneously  within  unconsolidated  porous 
media.  An  important  case  in  question  would  be  the  investigation  of 
water  evaporation  from  naturally  occurring  soils,  sands  and  gravels, 
since  this  phenomenon,  although  of  paramount  importance  in  civil 
engineering,  is  poorly  understood  at  present. 
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A. 1  ESSENTIAL  GEOMETRY  OF  THE  OBLATE  SPHEROID 

An  oblate  spheroid  is  generated  when  an  ellipse  is  rotated  about 
its  minor  axis  (Figure  26) .  The  system  of  coordinates  convenient  to 
this  investigation  will  be  oblate  spheroidal  coordinates  [5,n,ct>],  which 
constitute  the  special  case  of  ellipsoidal  coordinates  in  which,  of  the 
three  principal  axes  of  the  general  ellipsoid,  the  two  largest  are  of 
equal  length.  If  a  Cartesian  frame  of  reference  [x,y,z]  is  chosen  to 
be  collinear  with  these  three  principal  axes  then  it  can  be  shown^ 
that : 

x  =  a  sinh£  sinri  ,  (Al) 

y  =  a  cosh£  cosri  coscj)  ,  (A2) 

z  =  a  cosh£  cosn  sincj)  ,  (A3) 

where:  £  ^  0  ;  — tt /2  ^tt/2  ;  0  <  cj>  <  2tt  . 

In  these  equations  a  denotes  the  distance  between  the  trajectory  of  the 
focus,  F (0 ,a) ,  and  the  geometric  centre,  o.  The  coordinate  surfaces  of 
constant  £  comprise  a  family  of  confocal  oblate  spheroids;  those  of 

constant  n  comprise  a  family  of  confocal  hyperboloids  of  revolution. 

The  third  orthogonal  coordinate,  cj) ,  denotes  the  azimuthal  angle 
measured  out  of  the  plane  of  the  paper  of  Figure  26.  Typical  unit 
vectors  i  and  i  are  depicted  in  this  figure;  the  unit  vector  i,  is 

— £  — n  * 

orthogonal  to  the  paper. 
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FIGURE  26.  OBLATE  SPHEROIDAL  COORDINATES  IN  A  MERIDIAN  PLANE 
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An  element  of  distance  in  space,  dL,  may  be  represented  by: 


(dL)2  =  (dx)2  +  (dy)2  +  (dz)2  =  (hrd£)2  +  (h  dp)2  +  (h  d?)2  ,  (A4) 

^  n  9 

where  h^  ,  and  h^  denote  the  so-called  metrical  coefficients  (scale 

S  H  r 

factors)  which  follow  directly  from  Equations  (Al)  -  (A4) ,  viz. 


a  i^sinh  2g  +  sin11'  , 
a  cosh^  cosp  . 


(A5) 

(A6) 


The  particular  spheroidal  shape  of  interest  will  henceforth  be 
identified  by  c,  =  ;  for  this  geometry  it  nay  be  noted  that: 

length  of  seni-ninor  axis  =  a  sinhl^  ,  (A7) 

length  of  seni-najor  axis  =  a  coshl^  .  (A8) 

When  designating  specific  spheroidal  shapes  it  is  customary  to  do  so 
by  means  of  their  eccentricity,  defined  by: 


eccentricity  =  minor  axis /major  axis  -  e 

=  tanh^Q  0  <  e  <  1.0  .  (A9) 


This  definition  holds  also  for  prolate  spheroids,  and  will  be  used  as 
such  in  Appendix  B. 

For  future  convenience  the  following  functions  will  be  defined 


(A10) 

(All) 


f($)  =  1/sinh^  -  arccot(sinh-)  , 
g(r)  =  sinh^/cosh-^  -  arccot  (sinlu)  , 
hU)  «  2fU)  -  gU)  • 


(A12) 


. 


. 


' 
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A. 2  THE  PROPOSED  MODEL  FOR  CO-AXIALLY 
ORIENTATED  OBLATE  SPHEROIDS 

Attention  will  here  be  confined  to  an  homogeneous  swarm  of 
oblate  spheroids  of  identical  eccentricity,  e,  possessing  an  arbitrary 
size  distribution  and  orientated  such  that  the  axis  of  revolution  of 
each  particle  is  collinear  with  the  x-direction. 

The  proposed  model  for  this  system  (Figure  27)  is  closely 
related  to  that  for  spheres  (Figure  3) .  It  here  consists  of  a 
reference  spheroid  (£  =  £q) ,  an  annular  region  of  void  space  (bounded 
by  the  two  confocal  spheroidal  surfaces  £  =  £q  and  E,  =  E,^)  and  an 
exterior  region  of  homogeneous  porous  material  possessing  the  same 
macroscopically  anisotropic  characteristics  as  the  original  system. 

In  order  that  the  porosity  of  the  unit  cell  (comprising  the 
reference  spheroid  and  the  annular  region)  remains  equal  to  that  of 
the  original  system,  e,  it  is  necessary  that 

{  (4tt / 3)  (acosh^)2  (asinh£Q)  }/{  (4tt/3)  (acosh^)  2  (asinh^)  }  =  (1-e)  .  (A13) 

After  re-arrangement,  this  expression  yields  the  condition: 

sinh3£^  +  sinh£^  -  (cosh^Qsinh^Q)  /  (1-e)  =  0  .  (A14) 

It  should  be  noted  that  the  eccentricity  of  the  outer  surface  of  the 
unit  cell  will,  in  general,  be  different  from  that  of  the  reference 
spheroid . 

This  model  possesses  the  distinct  advantage  that  the  macro¬ 
scopically  homogeneous  and  anisotropic  characteristics  of  the  original 
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A. 2] 


system  of  spheroids  remain  unaffected  by  the  modelling  procedure,  i.e0 
the  modelled  system  is  macroscopically  indistinguishable  from  the 
original  system  when  viewed  from  the  outsider, 


. 
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FIGURE  27.  THE  PROPOSED  MODEL  FOR  AN  HOMOGENEOUS  SWARM  OF  CO-AXIALLY  ORIENTATED 
OBLATE  SPHEROIDS  (Cross-Section  Through  Centre  of  Reference  Spheroid) 
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A. 3  THE  DEFINING  EQUATIONS  AND  BOUNDARY  CONDITIONS 
A. 3.1  GENERAL  DISCUSSION 

The  macroscopic  form  of  Fick's  Law  which  describes  diffusion 
through  an  homogeneous  swarm  of  co-axially  orientated  oblate  spheroids 
may  be  recalled  from  Equations  (45)  -  (51)  of  the  main  text  as  being: 

q*  =  -D  A  Vc*  ,  (A15) 


with  the  tensor  A,  in  the  preferred  frame  of  reference  [x,y,z],  being 
represented  by 


A 


A  0  0 

x 

0  A  0 

y 

0  0  A 

y 


(A16) 


It  may  be  further  recalled  that  in  order  to  fully  describe  diffusion 

occurring  in  an  arbitrary  direction  it  suffices  to  possess  a  knowledge 

of  A^  and  A  ,  since  this  provides  all  the  necessary  information  to 

construct  the  tensor  A  appearing  in  Equations  (A15)  and  (A16)  above. 

Consequently,  the  principal  objective  of  this  investigation  must  be 

to  so  evaluate  A  (the  diffusivity  factor  in  the  x-direction)  and  A 

x  y 

(the  diffusivity  factor  in  the  y-z  plane) . 

Intrinsically,  the  modelled  system  possesses  the  same  macroscopic 

anisotropy  (characterized  by  A)  as  the  original  system  of  spheroids. 

However,  the  evaluation  of  A  and  A  for  the  modelled  system  implicitly 
»  x  y 

presupposes  a  knowledge  of  A  (which  itself  implies  a  knowledge  of  A 
and  A  ).  The  problem,  as  it  stands,  is  therefore  intractable. 


. 
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However,  this  difficulty  may  be  averted  for  diffusion  occurring 
specifically  in  a  characteristic  direction.  In  general,  the  mainstream 
flux,  Q* ,  within  the  original  system  will  not  be  collinear  with  its 
inducing  concentration  gradient  on  account  of  the  anisotropy  of  the 
system.  Only  for  diffusion  occurring  specifically  in  a  characteristic 
direction  will  Q*  be  collinear  with  its  inducing  concentration  gradient. 
In  other  words,  for  flow  in  a  characteristic  direction  the  macroscopic 
flowfield  will  remain  totally  unaware  of  any  anisotropy  in  the  system. 
Therefore,  the  previously  mentioned  difficulty  within  the  modelled 
system  may  be  averted  for  diffusion  occurring  specifically  in  a  charac¬ 
teristic  direction  by  treating  the  macroscopic  flowfield  within  the 
exterior  region  of  this  system  as  if  it,  too,  is  unaware  of  any 
anisotropy  of  the  system  as  a  whole+.  Thus,  Fick’s  Law  for  the 
exterior  region  ^  £  <  00 )  may  be  written  in  the  form: 

q*(5  ,n  ,<J>)  =  [qj£,q*,q|] 

=  [-D*(l/h  )(3c*/H)  ,-D*(l/hn)  Oc*/9n)  ,-D*(l/h  )  (3c*/3cf>)], 

(A17) 

where : 

D*  =  D*  for  diffusion  collinear  with  the  x-direction; 
x 

D*  =  D*  for  diffusion  collinear  with  the  y-z  plane. 

y 


t  It  will  be  recalled  from  Section  1.3.4  that  for  spheres  the 
exterior  flowfield  is  everywhere  collinear  with  the  mainstream  flowfield. 
Observations  made  in  Section  A. 4  to  follow  demonstrate  that  this  also 
holds  true  for  spheroids. 


„ 


.  . 


)&  !  .  ' 
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Fick’s  Law  for  the  isotropic  annular  region  (£q  <  £  <  is: 

q(5,n,<j>)  =  [q?,qn,q({)] 

=  [-D(l/h  )(3c/aO,-D(l/h  )  (9c/3n)  »~D  (1/h  )  (8 c/3 4> ) ]  .  (A18) 

c,  n  4> 

A. 3. 2  FUNDAMENTAL  DIFFERENTIAL  EQUATIONS 

By  direct  analogy  with  Equations  (7)  and  (8)  of  the  main  text 
(for  an  homogeneous  swarm  of  spheres)  the  describing  equations  for 
the  corresponding  case  of  co-axially  orientated  oblate  spheroids  will 
be : 


annular  region: 

o 

ii 

o 

CM 

> 

?0  "  5  "  Cl  . 

(A19) 

exterior  region: 

V2c*  =  0 

?!  <  C  <  “ 

(A20) 

9  J.  l 

The  Laplacian  operator,  V  ,  here  assumes  the  form  : 

V2  •=  (9/9sinq)  (cos2n(9/9sinri})  +  (9 /9sinh£) (cosh2?{9 /9sinh£ }) 

+  (h2/h2)  (92/3c|>2)  .  (A21) 

5  4> 

A. 3. 3  STIPULATED  BOUNDARY  CONDITIONS 

Again,  by  direct  analogy  with  Equations  (12 ) — (15 )  for  spheres, 
the  following  boundary  conditions  respectively  must  be  stipulated 
for  co-axially  orientated  oblate  spheroids,  viz. 

q5(5;.n.*)  =  0  >  (A22) 

c(?^,n,4>)  =  c*(£*,n,4>)  ,  (A23) 

q^(?1,n,(J>)  =  q|(£^»n»4>)  »  (A2A) 


Limit  q*(£,n,<J>)  =  Q*  > 

£  oo  — 


(A25) 


‘ 


where : 


q*  =  [q*,0,0]  £  or  diffusion  in  the  characteristic  x-direction; 

Q*  =  [0,Q*,0]  for  diffusion  in  the  characteristic  y-direction; 

Q*  =  [0,0,Q*]  for  diffusion  in  the  characteristic  z-direction, 

these  latter  three  expressions  referring  to  the  preferred  Cartesian 

frame  of  reference  [x,y,z]. 
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A. 4  SOLUTION  OF  THE  DEFINING  EQUATIONS 

A.  4.1  CASE  1:  DIFFUSION  COLLINEAR  WITH  THE  AXES  OF  REVOLUTION  OF  THE 
SPHEROIDS 

The  flowfield  in  this  case  will  be  collinear  with  the  charac¬ 
teristic  x-direction  and  will  therefore  exhibit  axial  symmetry;  this 
implies  that  q*  =  0  and  =  0  in  Equations  (A17)  and  (A18)  respectively. 
Consequently,  the  <j>  variable  may  hereafter  be  suppressed. 

Sufficiently  general  solutions  of  Laplace's  Equations  (A19)  and 
(A20)  for  diffusion  in  the  x-direction  can  be  extracted  from  fundamental 
series  solutions  developed  by  Lamb^,  viz. 

c(£,n)  =  [A^+A^arccot (sinh£) ]  +  [A^+A^f (?) ]sinh£sinn  <5  <  (A26) 

C*(£,n)  =  [A*+A*arccot(sinh£)]  +  [A*+A*f  (?) ]sinh£sinn  ^  <  £  <  °°  >  (A27) 

where  A.  and  A^  (i  =  1,2, 3, 4)  denote  arbitrary  constants.  The  parti¬ 
cular  solutions  which  satisfy  the  boundary  conditions  stipulated  in 
Equations  (A22)-(A25)  transpire  to  be: 

c(£,n)  =  c*(5,0)  -  A(aQ*/D*)  [  1— f  (5)/g(5Q)  ]sinh£sinri  ,  (A28) 

c*(£,n)  =  c*(£,0)  -  (aQ*/D*)[l-A*f (5) ]sinh?sinn  ,  (A29) 

X 

in  which  A  and  A*  are  defined  by 

Xxg(50){f(5l)"g(5l)} 

A  =  — - 


f(C1){ga0)-g(C1)}  -  Xxg C51Hg(50)-f (5^} 


(A30) 


a ik)  ,  I  >  ;  *  0  <W  I.’f  1  -  •  "K:  .  > 
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A*  = 


{gU0)-gttn)}  -  xj gcg-fa,)} 


0- 


f(e1){gC50)-g(51)}  -  xxga1){ga0)-f(?1)} 


(A31) 


the  functions  f(£)  and  g(£)  being  defined  by  Equations  (A10)  and  (All) 
respectively . 

Since  the  modelled  system  is  to  be  quantitatively  representative 
of  the  original  system  then  it  is  necessary  that  within  the  unit  cell 
the  average  flux  in  the  mainstream  direction  be  equal  to  that  of  the 
mainstream  itself.  This  necessitates  that 


{  1/it  (acoshC ,  )  2}j  q  (£,0)  2Tr(acosh£)  d(acosh^)  =  Q*  .  (A32) 

1  F  n 

This  equation  is  analagous  to  Equation  (22)  of  the  main  text  for 
spheres.  The  flux  component  q^(£,0)  appearing  in  Equation  (A32)  may 
be  evaluated  by  means  of  Equations  (A18)  and  (A26) ,  thus: 

q  (?,0)  =  { -D (1/h  )(9c/3n)>  =  (AQ*/X  ){l-f (5) /g (5n) >  •  (A33) 

n  n  @n=0  x  u 

The  substitution  of  this  expression  into  Equation  (A32)  can  be  shown  to 
yield : 

/  sinhCU-f  (S)/g(S0)  }d(sinh?)  =  (X^/ 2A)  cosh2^  ^  .  (A34) 

After  integration  and  manipulation  this  equation  yields: 

1  -  gtep/gUo)  -  ^X/A  •  (A35) 

Following  substitution  for  A  from  Equation  (A30)  the  ultimate  expres¬ 
sion  for  the  diffusivity  factor  in  the  characteristic  x  direction 


can  be  obtained,  viz. 
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AX  =  {g^o)-g(?i)}/{ga0)-f(^i)}  • 


(A36) 


The  parameters  £  and  £  are  defined  by  Equations  (A9)  and  (A14) 
respectively . 

The  substitution  of  the  above  expression  for  A  into  Equation 


x 


(A31)  serves  to  illustrate  that,  in  fact,  A*  =  0;  this  implies  that  the 
flowfield  within  the  exterior  region  of  the  modelled  system  is 
everywhere  totally  undisturbed  by  the  modelling  procedure  (compare 
with  the  like  conclusions  of  Section  1.3.4  for  spheres). 

A.  4.2  CASE  2:  DIFFUSION  ORTHOGONAL  TO  THE  AXES  OF  REVOLUTION  OF  THE 
SPHEROIDS 

Although  the  following  analysis  relates  specifically  to  diffusion 
occurring  in  the  characteristic  z-direction  (for  which  the  calculations 
are  simplest)  the  final  results  will  be  valid  for  diffusion  occurring 
in  any  direction  collinear  with  the  y-z  plane  (from  considerations  of 
symmetry) .  As  already  noted  in  Equation  (50)  of  the  main  text  this 
implies,  in  particular,  that: 


(A3  7) 


However,  in  such  cases  the  flowfield  will  no  longer  be  axi-symmetric 
but  will  be  ^-dependent.  Sufficiently  general  solutions  of  Laplace's 
Equations  (A19)  and  (A20)  for  diffusion  in  the  z-direction  can  be 
shown ^  to  be: 

c(5  ,n»4>)  =  [B1+B2arccot(sinhC)]  +  [B3+B4g(£)  ]cosh£cosnsin<}> 


(A38) 


,  S.  ,  .  ■  .  . 


I 


> 
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=  [B*+B*arccot  (sinh£)  ]  +  [B*+B^g  (£)  ]  cosh£  cos  rising 

Kl  <£<*>,  (A39) 

where  B_^  and  B*  (i  =  1,2, 3, 4)  denote  arbitrary  constants.  The  parti¬ 
cular  solutions  which  satisfy  the  boundary  conditions  stipulated  in 
Equations  (A22)-(A25)  transpire  to  be: 


c(£,r),4>)  -  c*(£,n,0)  -  B  (aQ*/D*)  [l-g  (£)  /h  (£q)  ]cosh£cosnsin<j)  ,  (A40) 


C,!c  (?  » ri ,  4> )  -  c*(^,n,0)  -  (aQ*/D*)  [l-B*g  (£)  ]cosh£cosr|sin<}>  ,  (A41) 


in  which  B  and  B*  are  defined  by 


Ayh^0){gai)"h(Cl)} 

g(?1){h(C0)-h(51) }  -  X  h(51){h(50)-g(C1)} 

(h(?0)-h(ei)}  -  Xy{h(50)-g(C1)} 

ga1){ha0)-h(c1)>  -  x  ha1){hu0)-ga1)> 


(A42) 


(A43) 


the  functions  g(£)  and  h(£)  being  defined  by  Equations  (All)  and  (A12) 
respectively . 

Again,  since  the  modelled  system  is  to  be  quantitatively 
representative  of  the  original  system  then  it  is  necessary  that  within 
the  unit  cell  the  average  flux  in  the  mainstream  direction  be  equal  to 
that  of  the  mainstream  itself.  This  here  necessitates  that 

h 

{1/tt  (asinh^)  (acosh^) }/  q^(£,n,0)  dA  =  Q*  ,  (A44) 


where  dA  denotes  an  annular  element  of  area  in  the  plane  cf>  =  0 ,  viz . 


dA  =  TTa2(sinh2£  +  cosh2£)d£  . 


(A45) 


The  flux  component  q^(C»n»0)  appearing  in  this  equation  may  be 
evaluated  by  means  of  Equations  (A18)  and  (A38) ,  thus : 


q  (?,n,0)  =  {-D(l/h  )(3c/8(J>)} 

^  ^  @<p=0 


(BQ*/Ay){l-g(C)/ha0)}  . 


(A46) 


The  substitution  of  Equations  (A46)  and  (A45)  into  Equation  (A44)  can 
be  shown  to  yield: 


/  (cosh2£  +  sinh20{l-g(C) /h(?n)  }dg  =  (A  /B)sinhC ,  cosh?,  .  (A47) 

r  u  y  I  i 


After  integration  and  manipulation  this  equation  reduces  to: 


1  -  MCP/MCq)  =  X  / B  .  (A48) 

Following  substitution  for  B  from  Equation  (A42)  the  ultimate  expres¬ 
sion  for  the  diffusivity  factor  in  any  direction  collinear  with  the 
y-z  plane  can  be  obtained,  viz. 


Ay  =  {h(^0)-h(C1)}/{h(C0)-g(51)>  •  (A49) 

The  parameters  and  E,  ^  are  again  defined  by  Equations  (A9)  and  (A14) 
respectively . 

The  substitution  of  the  above  expression  for  Ay  into  Equation 
(A43)  shows  that,  in  fact,  B*  =  0 ;  this  implies  (once  again)  that  the 
exterior  flowfield  within  the  modelled  system  is  everywhere  totally 
undisturbed  by  the  modelling  procedure  (compare  again  with  the  like 
conclusions  reached  in  Section  1.3.4  for  spheres). 
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B.l  ESSENTIAL  GEOMETRY  OF  THE  PROLATE  SPHEROID 


A  prolate  spheroid  is  generated  when  an  ellipse  is  rotated 
about  its  major  axis.  The  system  of  coordinates  appropriate  to  this 
investigation  will  be  prolate  spheroidal  coordinates  [£,n»<f>]*  This 
coordinate  system  constitutes  the  special  case  of  ellipsoidal 
coordinates  in  which,  of  the  three  principal  axes  of  the  general 
ellipsoid,  the  two  smallest  are  of  equal  length.  If  a  Cartesian 
frame  of  reference  [x,y,z]  is  chosen  to  be  collinear  with  these  three 
principal  axes,  then  it  can  be  shown^  that: 


x  =  a  coshC  sinn  , 
y  =  a  sinh£  cosn  cos<})  , 
z  =  a  sinhC  cosn  sin<{>  , 


(Bl) 

(B2) 

(B3) 


where:  E,  -  0  ; 


-tt/2  <  n  ^  tt  /  2  ; 


0  <  <J)  <  2tt  . 


,,  a  |u 


■ 
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B.2  THE  SOLUTION  FOR  PROLATE  SPHEROIDS 


Attention  will  here  by  confined  to  an  homogeneous  swarm  of 

prolate  spheroids  of  identical  eccentricity,  e,  possessing  an  arbitrary 

size  distribution  and  orientated  such  that  the  axis  of  revolution  of 

each  particle  is  collinear  with  the  x-direction. 

However,  rather  than  solving  this  problem  anew,  it  is  possible 

to  proceed  more  directly  by  making  avail  of  the  preceding  results 

developed  for  oblate  spheroids  together  with  a  standard  mathematical 

42 

transformation  .  Thus,  as  may  readily  be  verified,  if  a  is  replaced 
by  -ia  and  sinh£  by  icosh£  in  the  expressions  developed  in  Appendix  A 
for  oblate  spheroids  then  the  solutions  of  the  corresponding  problems 
involving  prolate  spheroids  will  be  obtained.  Execution  of  this 
transformation  on  Equations  (A36) ,  (A49)  and  (A10)-(A12)  ultimately 
yields  the  following  expressions: 


=  <g^o)'G(51)}/{G(50)'F(51)}  ’ 


(B4) 


Xy  =  {H(50)-H«1)}/{H(C0)-G(51)}  , 


(B5) 


in  which  the  functions  F(£),  G(£)  and  H(£)  are  defined  by 


F(£)  =  1/coshC  -  areacoth(coshC)  , 


(B6) 


G(£)  =  cosh£/sinh2£  -  areacoth (cosh?)  , 


(B7) 


H(5)  =  2F(5)  -  G(€)  . 


(B8) 


. 


I  '  - 


-  .  o  -:r  ,x 

•• 


E. 2~ 


- .a  a  parameter  i  =  5.2 

Equation  (A9j,  viz. 


However,  tie  parauaaar 

cosh r  E  ^  -  cc 

This  equation  ensures 
spheroids  is  equal  to 
Equation  (A1-)  for  obi 
derived  frcn  this  equa 
specified  natheratlcal 


sIM 


■ 
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C.l  REVIEW  OF  THE  DEFINING  EQUATIONS  AND  BOUNDARY  CONDITIONS 

C.1.1  FUNDAMENTAL  DIFFERENTIAL  EQUATIONS 

The  differential  equations  describing  steady,  incompressible 
creeping  flow  through  the  modelled  system  (Figure  16)  have  been 
presented  and  developed  in  Section  II. 4  of  the  main  text.  Solutions 
are  there  sought  of  Equations  (94)  and  (95)  within  the  annular  and 
exterior  regions  respectively,  viz. 

annular  region:  E2(E2i|j)  =0  R  <  r  <  s  ,  (Cl) 

exterior  region:  -(1/k)E2^*  +  E2(E2^*)  =0  S  <  r  <  00  .  (C2) 

In  these  equations  and  ip *  denote  the  streamf unctions  in  the  annular 

and  exterior  regions  respectively,  and  E2  denotes  the  Spherical 

23  A  2 

Harmonic  Operator  ’  ,  defined  by: 

E2  e  (92/9r2)  +  (1/ r2) (92/902)  -  (cot0/r2) (9/96)  .  (C3) 

C.l. 2  STIPULATED  BOUNDARY  CONDITIONS 

The  boundary  conditions  which  must  be  satisfied  have  been 
stipulated  in  Equations  (97)-(103),  viz. 


ur(R  ,0)  =  0  , 

(C4) 

uq(R+,0)  =  0  , 

(C5) 

u  (S",0)  =  u*(S+,0)  , 

(C6) 

u0(s",0)  =  u*(S+,0)  , 

(C7) 

t(s",0)  =  T*(S+,0)  , 

(C8) 

P(S”,0)  =  P*(S+,0)  , 

(C9) 

•  I 


'  N 
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Limit  u*(r,0)  =  U  *  ,  (CIO) 

r  ->■  oo 

where  U*  =  [u*,0,0]  in  Cartesian  coordinates  [x,y,z].  The  velocity 
components  u^,  u^  ,  u£  and  u*  appearing  above  are  defined  in  terms  of 
their  respective  streamfunctions  by  Equations  (92)  and  (93),  viz. 


ur  (r  ,9)  = 

-(l/r2sin0) (3^/30)  , 

(Cll) 

ii 

/  N 

CD 

5-i 
v _ ' 

CD 

3 

(l/rsin0)  (difj/dr)  , 

(C12) 

uj(r,0)  = 

-  (l/r2sin0)  (3ip*/30)  , 

(C13) 

uq (r  ,0)  = 

(l/rsin0)  (dilj*/dr)  . 

(C14) 

For  the  present  case  of  axial  symmetry,  t  and  x*  denote  the  non-zero 

shear  stress  components,  x  and  x*  ,  of  their  respective  stress 

ru  r  u 

28 

tensors  ,  viz. 

x (r ,0 )  =  -y[r3 (uQ/r)/3r  +  (1/r) (3ur/30) ]  ,  (C15) 

x*(r,0)  =  -y[r3 (u*/r)/3r  +  (1/r) (3u*/30) ]  .  (C16) 

o  r 

The  pressures  p  and  p*  are  defined  respectively  by  the  Navier-Stokes 

Equation  (84)  and  the  Brinkman  Equation  (85).  These  equations  yield, 
28 

on  expansion  ,  the  following  expressions  for  the  radial  pressure 
gradients  3p/3r  and  3p*/3r: 

3p/3r  =  y [ (3 2u^/3r2)+ (2/r) (3u^/3r)+(l/r2) (32u^/302)+(cot0/r2) (3u^/30) 

-(2ur/r2)-(2/r2) (9u0/30)-(2uecotS/r2) ]  , 


(C17) 
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3p*/9r  =  u[(92u*/9r2)+(2/r)  (9u*/9r)+(l/r2) (92u£/902)+(cot0/r2) (9u*/90) 

- (2u*/r2)  -  (2/ r2) (9u*/ 9  0)- (2u*cot0 / r2)- (u*/k) ]  . 
r  0  0  r 

(C18) 


The  pressures  p  and  p*  may  be  obtained  by  integration  of  these  equations 
with  respect  to  r. 


. 
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C • 2  SOLUTION  OF  THE  DEFINING  EQUATIONS 
C.  2.1  TEE  GENERAL  SOLUTIONS 

ine  following  analysis  derives  in  detail  the  results  presented 
and  utilized  in  Sections  II. 5  and  II. 7  of  the  main  text. 

The  Exterior  Rea  ion 

& 

The  boundary  condition  implicit  in  Equation  (CIO)  may  be 
rewritten  in  component  form  as  follows,  viz. 


Limit  u*(r,S)  =  -U*cos6  , 

X  ->  30 

(C19) 

Limit  u?(r,0)  =  U*sinS  . 

X  ->  » 

(C20) 

These  equations,  in  connection  with  Equations  (C13)  and  (C14)  ,  imply 
that : 


Limit  v* (r ,9)  =  (l/2)U*r2sin28  .  (C21) 

00 

Observing  the  nature  of  this  limiting  solution  for  p*  a  general 
solution  of  Equation  (C 2 )  is  sought  in  the  separable  form: 

(r ,8)  =  [u1  (r)  +  u>2(r)  +  u^Cr)  +  w^(r)]sin20,  (C22) 

where  u^(r)  denotes  any  function  of  r  such  that  _^(r)sin20  constitutes 
a  particular  solution  of  Equation  (C 2) . 

It  transpires  that  a  general  solution  of  the  above  form  is: 

V*(Xf9)  =  (<U*/2) [E/x  +  Fx2  +  Ge”‘ (l+x)/x  +  He ‘ (1-x) /x]sin2S ,  (C23) 

where  E,F,G  and  H  denote  arbitrary  constants  and  x  represents  the 


. 


X 


. 
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normalized  radial  coordinate  defined  by: 

X  (r)  =  r//IT  .  (C24) 

From  a  closer  inspection  of  Equations  (C21)  and  (C23)  it  follows 
necessarily  that  H  =  0  and  F  =  1,  whence  Equation  (C23)  becomes: 

^*(x,e)  =  (kU*/2) [E/x  +  x2  +  Ge"X(l+x)/x]sin2e  .  (C25) 

The  Annular  Region 

A  general  solution  of  Equation  (Cl)  ,  possessing  the  same 

2  8 

separable  form  as  Equation  (C22)  ,  can  be  shown  to  be  : 

iKx.e)  =  (kU*/2)[A/x  +  BX  +  Cx2  +  DX4]  sin20,  (C26) 

where  A,B,C  and  D  denote  arbitrary  constants. 

The  Global  Stream  function  Distribution 
The  entire  flowfield  around  the  reference  sphere  is  completely 
defined  by  the  streamfunctions  presented  in  Equations  (C25)  and  (C26)  , 
viz . 

<Kx,6)  =  (kU*/2)[A/x  +  Bx  +  Cx2  +  DX4]sin2e  a  <  x  <  6  ,  (C27) 

4>*(x,e)  =  (kU*/2)[e/x  +  X2  +  Ge  x (1+x) /x]sin20  3  <  X  <  00  •  (C28) 
The  parameters  a  and  3  represent  the  particular  values  of  X  defined  by 

a  =  x(R)  =  R/'/k*  ,  (C29) 

3  =  x(s)  =  S/v4T  . 


(C30) 


. 
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Now  that  the  general  form  of  the  global  streamf unction  has  been 
ascertained  the  global  velocity,  shear  stress  and  pressure  distri¬ 
butions  may  be  determined  directly,  as  follows. 

The  Global  Velocity  distribution 

The  radial  and  tangential  components  of  the  prevailing  flow- 
field  follow  at  once  from  Equations  (C11)-(C14)  and  Equations  (C27)- 


(C28)  ,  thus : 

u  (x,9)  =  -U*[A/x3  +  B/x  +  C  +  Dx2]cos9  ,  (C31) 

uj(x,0)  =  “U*[l  +  E/x3  +  Ge  X (1+x) /x3]cos0  ,  (C32) 

u„ (x  » 0 )  =  U*[-A/ 2X3  +  B/2X  +  C  +  2DX2]sin9  ,  (C33) 

0 

u*(x,0)  =  U*[l  -  E/2X3  -  Ge"X(l+x+X2)/2x3]sin9  .  (C34) 

0 

The  Global  Shear  Stress  distribution 


From  Equations  (C15)-(C16)  and  Equations  (C31)-(C34)  above  it 
may  be  shown  that: 

t  (x  >0)  =  -  (yU*/ /i<)  [3A/x4  +  3Dx]sin9  ,  (C35) 

x*(x,0)  =  -(uU*/v/k)[3E/x4  +  Ge”X(6+6x+3x2+X3)/2x4]sin9.  (C36) 

The  Global  Pressure  distribution 

From  Equations  (C17)— (C18)  and  Equations  (C31)— (C34)  it  is 

possible  to  show  that: 


p  (X  ,0)  =  p*  (X  »7T/2)  -  (uU*//IO[B/x2  +  10DX]cos9,  (C37) 
P*(X>9)  =  P*(x»"/2)  -  (uU*/v/k)[e/2x2  -  xlcos9  .  (C38) 
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C.2. 2  EVALUATION  OF  TEE  ARBITRARY  CONSTANTS 

The  stipulated  boundary  conditions  presented  in  Equations  (C4)- 
(C9)  may  be  re-expressed  in  terms  of  the  normalized  parameters  a  and 
3,  as  follows: 


ur (ot+,0)  =  0  , 

(C39) 

uQ(a+,e)  =  0  , 

(C40) 

ur (3~,e)  =  u* (3+,e)  , 

(C41) 

ue  (3~  ,0)  =  u*(B +,0)  , 

(C42) 

t(3-,0)  =  t*  (3+,0)  , 

(C43) 

p(3~,0)  =  p* (3+,0)  . 

(C44) 

It  transpires  that  these  six  conditions  contain  sufficient  information 
to  permit  the  determination  of  the  six  unknown  arbitrary  constants 
A,B,C,D,E  and  G.  Thus,  by  applying  these  conditions  to  the  respective 
global  distributions  presented  in  Equations  (C31)-(C38)  the  following 


six  independent  equations  may  be  generated: 

A  +  Ba2  +  Ca3  +  Da5  =  0  ,  (C45) 
A  -  Ba2  -  2Ca3  -  4Da5  =  0  ,  (C46) 
A  +  B32  +  C33  +  D35  =  33  +  E  +  Ge"3(l+3)  ,  (C47) 
A  -  B32  -  2C03  -  4D35  =  -23s  +  E  +  Ge"6(l+3+32)  ,  (C48) 
6A  +  6D35  =  6E  +  Ge"6(6+63+332+33)  ,  (C49) 
2B  +  2OD03  =  -203  +  E  .  (C50) 


'  x 
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Equations  (C45)-ivC50)  constitute  six  simultaneous  equations  in  the 
six  unknowns  A,B,C,D,E  and  Gj  the  solution  of  these  equations  can  be 
shown  to  be: 

A  =  6a3  (-238-735-1534-1583+334cx2-33a3+333a2+32a3) /J(a,3)  , 

B  =  6a  (638+213  8+45B4+4533-534ot2-503a2-3a5-a5) /J(a,3)  , 

C  =  3(-836-2835-6034-6033+533a3-532a3+33a5+3a5)/J(a,3)  , 

D  =  3 (484+483-333a+332a-3a3-a3) /J(a,3)  , 

E  =  2(-439-2438-6037-6036+938a+4537a-1036a3+126B6a-3035a3 

+21685a+934a5-6034a3+27084a-483a6+983a5 

-6033a3+27033a-68ae-6a6)/J(a,3) 

G  =  6e3(2036-2785a+533a3-9033a+2a6)/J(a,3)  , 

where : 

J (a ,3)  =  (-436-2435-18034-18033+935a+4534a-1033a3+1803sa 

-3032a3+93a5“4a8+9a5)  . 

C.2. 3  TEE  RESISTANCE  OFFERED  BY  THE  REFERENCE  SPHERE 

The  normal  resistance  force,  F^,  and  the  tangential  resistance 

force,  F  ,  offered  by  the  reference  sphere  may  be  calculated  by  inte- 

grating  the  local  pressure  p(a,0)  and  the  local  shear  stress  x(a,0) 

2  8 

respectively  over  the  entire  surface  thereof  ,  thus: 

IT 

F  =  27TR2/  [p(a,0)cos0]  sin0  d0  ,  (C58) 

r  0 

TT 

F  =  2ttR 2f  [-T  (a,0)sin0]  sin0  d6  . 

0  0 


(C51) 

(C52) 

(C53) 

(C54) 

,  (C55) 
(C56) 


(C57) 


(C59) 


' 


' 
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After  substituting  for  p(a,0)  and  i(a,0)  from  Equations  (C37)  and  (C35) 
respectively,  it  is  found  that: 


Fr  =  67ryU*R(-2B-20Da3)/9a  , 

(C60) 

Fg  =  67TiiU*R(4A+4Da3) /3a3  . 

(C61) 

The  total  resistance  force,  F,  offered  by  the  reference  sphere  will 
therefore  be: 


F  =  Fr  +  F0  =  67ryU*R(12A-2Ba2-8Da5)/9a3  .  (C62) 

After  substituting  for  A,B  and  D  from  Equations  (C51) ,  (C52)  and 
(C54)  respectively  the  ultimate  expression  for  F  is  obtained,  viz: 

F  =  6ttuU*R  £(a,3)  >  (C63) 

where : 

£(a,3)  =  4(-636-2135-4534-4533+534a2+53V+3a5+a5)/J(a,3)  .  (C64) 

It  has  been  demonstrated  in  the  main  text  (Table  3)  that  as  e  -*  1.0, 

28 

5 (a, 3)  1.0  so  that  Equation  (C63)  then  approaches  Stokes  Law  for 

creeping  flow  past  a  single  isolated  sphere,  viz. 

FStokes  “  6^U*R  *  (C65) 

C.2. 4  THE  TOTAL  RESISTANCE  OFFERED  BY  AN  ARBITRARY  SWARM  OF  SPHERES 

A  randomly  textured  swarm  of  spheres  possessing  the  following 


size  distribution  will  be  considered: 


' 


C.2] 
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Total  number  of  spheres  :  N  , 

Number  of  spheres  having  radius  R±  :  N  (i  =  1,2,3,  . . . ,n) , 

Fraction  of  spheres  having  radius  R.  :  f.  =  N  /N, 

111 

Size  Ratio  :  y.  =  R./Rn, 
y  1  l  0 

where  R^  denotes  the  radius  of  any  chosen  species  of  sphere  present, 
for  example  the  radius  of  the  most  prolific. 

The  Original  System 

This  analysis  has  heretofore  considered  an  unbounded  swarm  of 
spheres.  This  implies  an  absence  of  all  wall  and  end  effects  so  that 
U*  will  be  uniformly  constant  throughout  and,  in  consequence, 

V2U*  =  0,  Under  these  circumstances  the  Brinkman  Equation  (82)  reduces 
to  the  familiar  Darcy  Equation  (81)  as  previously  noted  in  Section  II. 2 
of  the  main  text.  This  means  that: 

(u/k)U*  =  Ap*/Ax  ,  (C66) 

where  Ap*  denotes  the  pressure  differential  across  the  system  and  Ax 
its  overall  length.  Letting  A  denote  the  cross-sectional  area  of  the 
system  orthogonal  to  the  flowfield,  the  total  resistance  offered  by 
the  swarm  of  spheres  will  be: 

F  =  Ap*  A  .  (C67) 

swarm 

Substituting  for  Ap*  from  Equation  (C66)  then  yields: 

F  =  pU*AAx/<  =  yU*Va2/R2  #  (C68) 

swarm  u  u 

In  this  equation  V  denotes  the  total  volume  of  the  system  (AAx)  and  k 


.  ' 


■ 
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has  been  replaced  according  to  Equation  (C 29 )  ,  viz. 


K  =  Rq/«o  •  (C69) 

The  total  volume  of  the  system  may  be  expressed  in  terms  of  the  total 
particle  volume  and  the  porosity  as  follows: 


V  =  £(4tt/3)  r3n  /(1-e)  =  (4tt/3)RqN  Jjy^f  /(1-e)  .  (C70) 

i  i 

Substituting  this  expression  for  V  into  Equation  (C68)  then  yields 
the  following  expression  for  the  total  resistance  offered  by  the 
original  swarm  of  spheres,  viz. 

F  =  6TTpU*RnN{  2a^/9  (1-e)  }Ty?f .  .  (C71) 

swarm  0  0  “  l  l 

l 


The  Modelled  System 

According  to  Equation  (C63)  the  resistance  force,  F  ,  offered 
by  a  reference  sphere  of  radius  R^  will  be: 


F.  =  6ttpU*R.  £  (a.,  3.)  =  6ttpU*R  y  5  (a  ,3.)  , 


(C72) 


where : 


a .  =  R.  / , 

l  l 

(C73) 

6  .  =  S  .  //<  . 
l  i 

(C74) 

There  still  exists  a  constant  ratio  between  &  and  ai  according  to 
Equation  (116)  of  the  main  text,  since  each  sphere  within  the  swarm 
may  be  considered  to  be  the  reference  sphere  in  turn,  thus. 


C.2] 


140 


=  •  (C75) 

The  total  resistance  offered  by  the  swarm  of  N  spheres  according  to 
the  proposed  model  may  therefore  be  expressed  by: 

Fm0del  “  l  FiNi  ’  N  l  Fifi  •  <C76> 

1  1 

Combining  Equations  (C76)  and  (C72)  then  yields  the  ultimate 
expression : 


Fmodel  =  I  •  «”) 

i 

C.  2.5  QUANTITATIVE  CONSISTENCY  OF  TEE  MODELLED  SYSTEM 

Since  the  modelled  system  is  to  be  quantitatively  representative 
of  the  original  system  (as  regards  hydrodynamic  resistance)  then  the 
total  resistances  as  predicted  by  Equations  (C71)  and  (C77)  must  be 
identical,  i.e. 


F 


model 


F 

swarm 


(C78) 


Equating  these  expressions  finally  generates  the  identity: 


2a0  ^  yifi  “  I  yifi^ai,ei^  =  0  *  (c79) 

i  1  i 

The  parameters  a.  and  may  be  expressed  in  terms  of  y^,  cXq  and  e 
as  follows,  thus: 


a  =  Ri/v/iT  =  (Rq/ZO  (R./R0)  =  aQy.  , 


(C80) 


ei =  1/3  =  °oyi(1’e)  1/3 


(C81) 


.  J  ‘  O'J  3  i 


X 
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After  introducing  these  expressions  into  Equation  (C79)  it  becomes 
apparent  that  otQ  and  e  are  related  by  an  exceedingly  complicated 
function,  viz. 

2ao  £  y|fi  "  1  y±f±^(<CLQy±»a0y±^-^)~1^3)  =  o  .  (C82) 

i  i 

Provided  the  size  distribution  (f^,y^)  of  the  spheres  is  known  then 
otQ  may  be  extracted  from  Equation  (C82)  as  a  function  of  e  by  an 
iterative  procedure. 


C.2. 6  COMPARISON  OF  THE  TOTAL  RESISTANCE  WITH  THAT  PREDICTED  BY 
STOKES  LAW 

The  total  resistance  of  the  original  swarm  of  N  spheres,  assum¬ 
ing  each  sphere  to  be  hydrodynamically  independent  of  the  remainder, 
would  be  given  by  Stokes  Equation  (C65)  as: 


Stokes 


£67ryU*RJSL  =  6TTyU*RQN  . 


(C83) 


It  is  instructive  to  compare  this  expression  with  the  resistance  pre¬ 
dicted  by  the  presented  theory.  Thus,  from  Equations  (C71)  and  (C83) 
it  follows  that: 


W  =  F  ,  /F  =  {9(l-e)/2o^}Qy  f  .]/[[y?f .]  .  (C84) 

Stokes  swarm  0  J  i  i  “  i  i 

Hence,  for  a  specified  porosity  e  and  a  known  size  distribution 
(f.,y<),  a  knowledge  of  aQ  is  sufficient  to  evaluate  W. 


Monosized  Spheres 

For  a  system  composed  exclusively  of  monosized  spheres 
(n=l,  i  =  1,  f.  =1,  yi=l  and  aQ  E  a)  Equations  (C82)  and  (C84) 


^  -V 


■  ->*v 
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reduce  to : 

2a2  -  9 (1-e) £ (a ,a(l-e)-1 / 3)  =  0  , 

W  =  9(l-e)/2a2  =  l/£(a,a(l-e)-1/3)  . 

Therefore,  for  a  specified  porosity  e  the  corresponding  value  of 
may  be  obtained  from  Equation  (C85)  by  an  iterative  procedure; 

W  may  then  be  computed  directly  from  Equation  (C86) . 


(C85) 

(C86) 

a 
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APPENDIX  D 


INCOMPRESSIBLE  CREEPING  FLOW  THROUGH 
A  FRACTURED  POROUS  MEDIUM 
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D.l  DISCUSSION  OF  THE  SYSTEM 

Fractures  in  petroleum  reservoirs  often  play  an  important  role 

in  the  production  of  fluid  from  the  pay  formation^.  In  fact,  in  some 

instances  the  oil  or  gas  could  not  be  produced  economically  in  the 

absence  of  fractures  functioning  as  arteries  for  transporting  the 

fluids  from  the  reservoir  rock  formation  to  the  well-bore.  Indeed, 

deliberate  attempts  are  often  made  to  induce  fractures  and  to  widen 

39 

existing  fractures  during  oil-field  production 

The  idealized  system  depicted  in  Figure  28  may  be  considered 
representative  of  a  fractured  petroleum  reservoir.  This  system  com¬ 
prises  two  semi-infinite  extents  of  isotropic  porous  medium,  possess¬ 
ing  permeabilities  and  (which  may  or  may  not  be  equal) , 
separated  by  a  uniform  fracture  (channel)  of  width  2h.  The  flowfield 
will  here  be  presumed  to  be  collinear  with  the  fracture. 

The  analysis  which  follows  seeks  to  calculate  the  velocity  pro¬ 
file  prevailing  throughout  the  idealized  system  described  above  and, 
in  doing  so,  to  demonstrate  the  inadequacy  of  the  Darcy  Equation  for 
describing  fluid  flow  in  porous  media. 


. 
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FIGURE  28. 


AN  IDEALIZED  FRACTURE  WITHIN  A  PETROLEUM  RESERVOIR 


• 
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D . 2  THE  DEFINING  EQUATIONS  AND  BOUNDARY  CONDITIONS 

D.2.1  FUNDAMENTAL  DIFFERENTIAL  EQUATIONS 

The  Brinkman  Equation  (82)  will  be  employed  to  describe  the 
hydrodynamic  conditions  within  each  of  the  porous  regions,  whilst  the 
Navier-Stokes  Equation  (80)  will  be  acknowledged  within  the  separating 
channel  of  free  fluid  space,  thus: 

-(y/x^)u*  +  y(d2u*/dy2)  =  dp*/dx  h  <  y  <  00  ,  (Dl) 

y(d2u/dy2)  =  dp/dx  -h  <  y  <  h  ,  (D2) 

-(p/<2)u*  +  y(d2u*/dy2)  =  dp*/dx  <  y  <  -h  .  (D3) 

At  great  distances  from  the  fracture  the  velocity  profiles  will  become 
uniformly  constant,  that  is: 

Limit  u*(y)  =  U*  ,  (D4) 

y  ->  oo 

Limit  u*(y)  =  U*  .  (D5) 

y  ->  — oo 

Substituting  these  expressions  into  Equations  (Dl)  and  (D3)  yields: 

-(y/<1)U*  =  dpj/dx  ,  (D6) 

-(y/K2)U*  =  dp*/ dx  .  (D7) 


However,  from  considerations  of  equilibrium  it  is  necessary  that 

dp*/dx  =  dp/dx  =  dp*/dx  =  constant  .  (D8) 

r  1  ^ 


N 


5>V  »  • 
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In  view  of  these  equalities  the  substitution  of  Equations  (D6)  and 
(D7)  into  Equations  (Dl)-(D3)  ultimately  yields: 

— u J/ k ^  +  d2u*/dy2  =  -UJ/k^  h  <  y  <  00  ,  (D9) 

d2u/dy2  =  -U*/k  -h  <  y  <  h  ,  (DIO) 

-u^/k^  +  d2u*/dy2  =  -U*/^  -«>  <  y  <  -h  .  (Dll) 

D.2. 2  STIPULATED  BOUNDARY  CONDITIONS 

From  considerations  of  continuity  and  equilibrium  at  the 
extrema  of  the  channel  the  following  boundary  conditions  must  be 
stipulated : 


u*(h+)  =  u(h  )  , 

(D12) 

Tj(h+)  =  T(h")  , 

(D13) 

u(-h+)  =  u*(-h  )  , 

(D14) 

x(-h+)  =  T*(-h  )  . 

(D15) 

In  these  equations  t*,  t  and  t*  denote 

the  tangential  shear  stresses, 

defined  by: 

t*  =  -y(du*/dy)  , 

(D16) 

t  =  -y(du/dy)  , 

(D17) 

t*  =  -y(du*/dy)  . 


(D18) 


' 
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D . 3  SOLUTION  OF  THE  DEFINING  EQUATIONS 

The  solutions  of  Equations  (D9)  -  (Dll)  which  satisfy  the 
boundary  conditions  stipulated  in  Equations  (D12)-(D15)  transpire  to 
be: 


u*(n)  =  U*(l+AJe_n)  6  <  n  <  00  ,  (D19) 

u(n)  =  u*(B+Cn-n2/2)  -6  <  n  <  <5  ,  (D20) 

u*(n)  =  U*(l+A*eYn)  _oo  <  n  <  -6  ,  (D21) 

U*  =  U*/y2  .  CD 22) 

In  these  equations  x\  denotes  the  normalized  variable  defined  by: 

n  =  y//ic^  .  (D23) 

The  parameters  6  and  y  are  defined  by: 

6  =  h//ic^  ,  (D24) 

y  =  /k^/i<2  >  (D25) 

and  A*  B,  C  and  A*  possess  the  specific  forms: 

_L  ** 

A*  =  (l-Y2+26y+262Y2) /M(6 ,y)  >  (D26) 

B  =  (1+t+6+26y+<5y2+362y/2+362y2/2+63y2) /H(6  ,y)  ,  (D27) 

C  =  -(1-y2+6y-6y2)/M(6,y)  ,  (D28) 

A*  =  -ye6Y(l-Y2-26Y2-262Y2)/M(6,Y)  > 


(D29) 


. 
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where : 


M(6,y)  =  y(1+Y+26y)  .  (D30) 

Equations  (D19)  -  (D30) ,  being  of  a  rather  general  nature,  are 
capable  of  several  interesting  specializations,  as  follows. 

One  Medium  Impermeable 

Supposing  that  medium  2  is  impermeable  ( < ^  =  0,  y  =  00 )  then  the 
prevailing  flowfield  will  be  described  by: 

u*(n)  =  U*(l+A*e~n)  ,  (D31) 

u(n)  =  U*(B+Cn-n2/2)  ,  (D32) 

u*(n)  =  0  ,  (D33) 

where : 

A*  =  -e6(l-262)/(l+26)  ,  (D34) 

B  =  6  (1+36/2+62)  / (1+26)  ,  (D35) 

C  =  (1+6)/ (1+26)  .  (D36) 

From  these  equations  it  is  clear  that  u(-6)  =  0 ,  as  required  for  non¬ 
slip  fluid  flow  past  a  solid  surface. 

Both  Media  Possessing  an  Equal  Permeability 

Under  these  conditions  (<2  =  Kl*  Y  =  1)  the  flowfield  will  be 


described  by: 


.  ‘  •'  •  -  H  J 
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U*(n)  =  U*(l+Se6_n)  ,  (D37) 

u(n)  =  U*(l+6+62/2-q2/2)  ,  (D38) 

u*(n)  =  UjCl+Se6^)  .  (D39) 

The  flowfield  is  therefore  symmetric  about  n  =  0,  as  would  be  expected. 

Both  Media  in  Contact 

In  this  case  (h  =  0,  6  =  0)  the  channel  separating  the  two  media 
ceases  to  exist  and,  in  effect,  constitutes  a  discontiriuity  between 
them.  The  flowfield  will  then  be  described  by: 


u*(n)  =  U*(l+A*e~n)  ,  (D40) 

u*(n)  -  (U*/y2)(l+A*eYn)  ,  (D41) 

where : 

A*  =  (1-y)/y  .  ^D42) 

A*  =  (y-1)  .  (D43) 


Furthermore,  when  the  permeabilities  of  the  two  media  become  equal  the 
flowfield  will  become  uniformly  constant.  Under  these  conditions 
(y  =  1)  Equations  (D40)  and  (D41)  reduce  to  the  necessary  forms,  viz. 

u*(n)  =  u*(n)  =  UJ  . 


(D44) 
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D.4  COMPARISON  OF  THE  PRESENTED  SOLUTION  WITH  THE  APPROXIMATE 
REPRESENTATION  OBTAINED  USING  THE  DARCY  EQUATION 

When  estimating  liquid  flow  through  fractured  porous  media  it 
44 

has  been  customary  in  the  past  to  adopt  the  Darcy  Equation  (81)  within 

28 

the  porous  regions  and  the  Poiseuille  Equation  within  the  separating 
channel.  Thus,  when  both  media  possess  equal  permeabilities  (which  is 
usually  the  case  in  practice)  the  flowfield  would  be  approximated  by: 


u*(n)  =  u* 

6  <  r)  <  00  , 

(D45) 

u(n)  =  u*(62/2-n2/2) 

-6  <  n  <  <5  , 

(D46) 

u*(n)  =  u* 

-CO  <  f|  <  -6  . 

(D47) 

These  profiles  are  depicted  schematically  in  Figure  29.  It  is 
immediately  apparent  that  the  above  equations  do  not  satisfy  the 
continuity  boundary  conditions  implicit  in  Equations  (D12)  and  (D14)  , 
viz . 


uj(6)  =  u(6)  , 

(D48) 

u(-6)  =  u*(-S)  . 

(D49) 

Moreover,  they  are  also  unable  to  satisfy  the 

implicit  in  Equations  (D13)  and  (D15) ,  viz. 

equilibrium  conditions 

x*(6)  =  x(<5)  , 

(D50) 

t(-6)  =  t*(-6)  . 


(D51) 


■ 


■ 
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RIGOROUS  SOLUTION 
APPROXIMATE  SOLUTION 


FIGURE  29.  THE  APPROXIMATE  AND  RIGOROUS  VELOCITY  PROFILES 
FOR  THE  FRACTURED  SYSTEM  (Schematic) 
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In  order  to  obtain  an  estimate  of  the  relative  error,  A, 
incurred  by  adopting  this  approximate  representation  of  the  flowfield 
it  will  be  instructive  to  compare  the  mean  velocity,  u,  within  the 
channel  as  predicted  rigorously  from  Equation  (D38)  with  that  pre¬ 
dicted  approximately  from  Equation  (D46)  ,  viz. 


A  =  (u  -  u  )/u  , 

approx 


(D52) 


where : 


+6 


u  =  (1/26)/  U*(l+6+62/2-n2/2)  dn  =  U*(l+6+62/3)  , 


-6 


+6 


u 


approx 


=  (1/26)/  U*(62/2-n2/2)  dn  =  U*(62/3)  . 

-6  1 


(D53) 

(D54) 


Substituting  Equations  (D53)  and  (D54)  into  Equation  (D52)  yields 


A  =  (l+6)/(l+6+62/3)  . 


(D55) 


This  equation  predicts  that  for  6  <  300  the  relative  error  A  will  be 
greater  than  1%.  The  immediate  inference  is  that  the  incurred  error 
will  be  of  practical  significance  for  sufficiently  thin  fractures. 
However,  the  thickness  (2h)  of  a  typical  petroleum  reservoir  fracture 
would  be  of  the  order  of  1  mm.,  with  6  being  in  the  range  500-5000. 

In  practice,  therefore,  the  quantitative  error  incurred  by  employing 
the  Darcy  Equation  (81)  rather  than  the  Brinkman  Equation  (82)  will 
usually  be  negligible. 

This  analysis  has  indeed  served  to  demonstrate  that  the  Brink- 
man  Equation  permits  a  mathematically  consistent  solution  to  a  realistic 
problem  involving  adjacent  regions  of  porous  medium  and  fluid  space;  a 
problem,  in  fact,  which  could  not  otherwise  be  treated  rigorously. 


! 
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APPENDIX  E 


INCOMPRESSIBLE  CREEPING  FLOW  THROUGH 
AN  ISOTROPIC  POROUS  MEDIUM 


CONTAINING  A  SPHERICAL  CAVITY 
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E.l  DISCUSSION  OF  THE  SYSTEM 

In  many  naturally  occurring  rock  formations  there  occur  distinct 
cavities  of  various  shapes  and  sizes.  Having  particular  significance  in 
Canada  are  those  encountered  in  oilfields  formed  within  dolomite  reefs. 
An  understanding  of  liquid  flow  through  such  two-porosity  systems  is 
therefore  desirable.  In  order  to  gain  physical  insight  into  this 
problem  it  will  be  instructive  to  consider  incompressible  creeping  flow 
through  the  idealized  system  consisting  of  an  unbounded,  isotropic 
porous  medium  which  contains  an  isolated  spherical  cavity.  As  yet, 
there  does  not  appear  to  be  available  a  solution  to  this  important 
problem. 

This  analysis  will  also  serve  to  demonstrate  the  necessity  and 
effectiveness  of  the  Brinkman  Equation  (82)  in  solving  another, 
hitherto  intractable,  two-region  problem. 


156 


E . 2  THE  DEFINING  EQUATIONS  AND  BOUNDARY  CONDITIONS 

E.2.1  FUNDAMENTAL  DIFFERENTIAL  EQUATIONS 

The  Navier-Stokes  Equation  (80)  will  describe  the  flowfield 
within  the  cavity  (of  radius  S) ;  this  equation  may  be  expressed  in 
terms  of  the  streamf unction  \p  as  per  Equation  (94),  viz. 

E2(E2^)  =0  0  <  r  <  S  .  (El) 

§ 

The  Brinkman  Equation  (82)  will  be  acknowledged  to  describe 
the  flowfield  within  the  isotropic  porous  medium  (of  permeability  k) ; 
in  terms  of  the  streamf unction  ip*  this  equation  assumes  the  following 
form  (Equation  (95)): 

-(1/k)E2i|j*  +  E2(E2^*)  =0  S  <  r  <  00  .  (E2) 

E.2.2  STIPULATED  BOUNDARY  CONDITIONS 

The  following  boundary  conditions  must  be  stipulated  for  this 


physically  realistic  problem: 

p(S~,0)  =  p*(S+,6)  ,  (E3) 

t(s",0)  =  x*(S+,0)  ,  (E4) 

u  (S',0)  =  u*(S+,0)  ,  (E5) 

uQ  (S~  ,0)  =  u*(S+,0)  ,  (E6) 

Limit  _u*(r,0)  =  U_*  .  (E7) 

oo 


The  arguments  on  which  the  above  boundary  conditions  are  based  are 
presented  in  Section  II. 4. 2  of  the  main  text.  In  addition,  a 
restriction  demanding  finiteness  of  the  flowfield  must  be  imposed, 
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in  particular  that 


ur(O,0)  must  be  finite  . 


(E8) 
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E.3  SOLUTION  OF  THE  DEFINING  EQUATIONS 

E.3.1  THE  GLOBAL  STREAMFUNCTI ON  DISTRIBUTION 

The  solutions  of  Equations  (El)  and  (E2)  which  satisfy  the 

boundary  conditions  implicit  in  Equations  (E3)-(E8)  can  be  shown  to 
be : 


^(x,9) 

-  (<U*/2)[Cxz  +  Dx4]sin^0 

0  <  X  <  B  ,  (E9) 

(x»0) 

=  (<U*/2>[e/X  +  X2  +  Ge"' (l+X)/x]sin20 

B  <  x  <  00  ,  (E10) 

where : 

X  =  t/Jk  , 

(Ell) 

8  =  S  /  Jk  , 

(E12) 

and: 

C  =  (663+2162+456+45)/(63+6B2+456+45)  , 

(E13) 

D  =  -3(B+1)/(B3+6B2+45B+45)  , 

(E14) 

E  =  2B3  (B3+6B2+15B+15)/(B3+6B2+45B+45)  , 

(E15) 

G  =  -30B3e?/(B3+6B2+45B+45)  . 

(E16) 

When  the  cavity  ceases  to  exist  (S  =  0 ,  6=0),  Equation  (E10) 
reduces  to  the  limiting  form  \p*  =  (l/2)U*r2sin20 .  This  is  the  familiar 
streamfunction  representation*-  of  an  undisturbed  flowfield,  U* ,  as 


would  be  expected. 
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E.3.  2  GENERAL  APPEARANCE  OF  THE  STREAMLINES 

Typical  streamlines,  as  computed  from  Equations  (E9)-(E16) 
above,  are  depicted  in  Figure  30  for  the  representative  case  3  =  1000. 
The  most  important  inference  to  be  drawn  therefrom  is  that  the  dis¬ 
turbance  created  by  the  cavity  on  the  mainstream  flow  is  effectively 
localized  to  a  region  concentric  with  the  cavity  and  possessing 
thrice  its  radius. 

It  is  evident,  therefore,  that  the  Brinkman  Equation  (82)  per¬ 
mits  the  solution  of  a  realistic  hydrodynamic  problem  which  has  not 
hitherto  been  solved. 

E.3. 3  THE  RRE SEN TED  SOLUTION  AS  A  LIMITING  SPECIALIZATION 

The  geometric  model  proposed  in  Section  II. 3  of  the  main  text 
for  an  homogeneous  swarm  of  spheres  (Figure  16)  is  fundamentally 
related  in  structure  to  the  spherical  cavity  system  here  investigated. 
The  modelled  system  may,  in  fact,  be  constructed  from  this  cavity 
system  by  suspending  a  sphere  (of  radius  R)  concentrically  within  the 
cavity . 

As  might  be  expected,  the  presented  solution  for  the  spherical 
cavity  system  agrees  identically  with  the  limiting  solution  for  the 
modelled  system  (Equations  (104)-(115))  for  the  limiting  specializa¬ 
tion  in  which  the  radius,  R,  of  the  reference  sphere  becomes  equal 


to  zero. 
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